# Data-Flow-Based Normalization Generation Algorithm of R1CS for Zero-
Knowledge Proof

Chenhao Shi, Hao Chen, Ruibang Liu, Guoqiang Li † Corresponding author.
Shanghai Jiao Tong University, Shanghai 200240, China
{undefeated, Dennis.Chen, 628628<EMAIL_ADDRESS>
###### Abstract

The communities of blockchains and distributed ledgers have been stirred up by
the introduction of zero-knowledge proofs (ZKPs). Originally designed to solve
privacy issues, ZKPs have now evolved into an effective remedy for scalability
concerns and are applied in Zcash (internet money like Bitcoin). To enable
ZKPs, Rank-1 Constraint Systems (R1CS) offer a verifier for bi-linear
equations. To accurately and efficiently represent R1CS, several language
tools like Circom, Noir, and Snarky have been proposed to automate the
compilation of advanced programs into R1CS. However, due to the flexible
nature of R1CS representation, there can be significant differences in the
compiled R1CS forms generated from circuit language programs with the same
underlying semantics. To address this issue, this paper uses a data-flow-based
R1CS paradigm algorithm, which produces a standardized format for different
R1CS instances with identical semantics. By using the normalized R1CS format
circuits, the complexity of circuits’ verification can be reduced. In
addition, this paper presents an R1CS normalization algorithm benchmark, and
our experimental evaluation demonstrates the effectiveness and correctness of
our methods.

###### Index Terms:

Zero-knowledge proof; Rank-1 constraint systems; Data flow graph; ZKP
programming; Normalization

##  I Introduction

_Zero-knowledge proofs (ZKPs)_ are increasingly recognized for their
importance in modern cryptography [1], as one and more cryptographic
communities seek to address some of the blockchain’s most significant
challenges: privacy and scalability. It is also the essential technique in
Zcash [2, 3]. From both user’s and developer’s perspectives, the heightened
emphasis on information privacy and security has led to a greater appreciation
for the privacy advantages offered by zero-knowledge proofs. As decentralized
finance (DeFi) usage grows, zero-knowledge applications that provide
scalability and privacy advantages will have more opportunities to increase
industry-wide adoption. However, not all computational problems can be
directly addressed using zero-knowledge proofs. Instead, we must transform the
issue into the correct form of computation. The _rank-1 constraint system
(R1CS)_ describes the execution of statements written in high-level
programming languages and is used by many ZKP applications, but there is no
standard way of representing them [4]. Circom is a novel domain-specific
language for transforming computational problems into R1CS format circuits.[5]
In the specific process of a first-order zero-knowledge proof, we first
convert the problem into a computational problem in Circom, then into R1CS
format circuits.

Due to the flexible nature of R1CS representation and variation of program
organizations and compiler optimization levels, there can be significant
differences in the compiled R1CS forms generated from circuit language
programs with the same underlying semantics, which leads to difficulties in
further ZKP program analysis and verification.

This paper proposes a data-flow-based algorithm for generating normalization
of R1CS, enabling the conversion of different R1CS constraints into a normal
form, facilitating the determination of equivalence and correctness. To
achieve this, the algorithm starts by transforming an R1CS into a data flow
graph structure resembling an expression tree. It then segments and abstracts
the data flow graph, eliminating differences between equivalent R1CS
constraints that may arise from the generation process. Finally, sorting rules
are proposed to sort the constraints and variables within R1CS, ultimately
resulting in a unique normal form for equivalent R1CS.

Moreover, we classify and summarize the reasons and characteristics of the
different equivalent R1CS generated, based on the constraint generation logic
of mainstream compilers and the expressiveness of R1CS. In addition, based on
the identified reasons for producing equivalent R1CS, we create a relatively
complete benchmark. Our proposed algorithm, which can pass all test cases in
the benchmark, demonstrates that equivalent R1CS can be converted into a
unique and identical canonical form under various circumstances.

This work contributes to R1CS optimization by providing a novel algorithm for
generating canonical forms of equivalent R1CS constraints. Our algorithm can
eliminate unnecessary redundancy and normalize representation, thus it can
improve existing methods and facilitates the analysis of equivalence and
correctness. Furthermore, the effectiveness and practicality of the proposed
algorithm are demonstrated through our comprehensive benchmark.

Related Work Eli et al. design, implement, and evaluate a zero-knowledge
succinct non-interactive argument (SNARG) R1CS [6]. Historically, research
investigating the factors associated with R1CS has focused on satisfiability.
In paper[7], where the prover incurs finite field operations to prove the
satisfiability of an n-sized R1CS instance [7]. Alexander et al. introduce
Brakedown, the first built system that provides linear-time SNARKs for NP [8].

Considering Circom and R1CS format circuits as two languages before and after
compilation, research on the generation of the R1CS paradigm is more akin to
research on semantic consistency in the compilation. Currently, the patent
applications and research papers propose ideas and solutions for generating
compilation paradigms in other languages, mainly exploring data flow [9],
syntax tree [10], or semantic mapping [11] aspects. These studies offer
crucial insights into the fundamental information semantically identical
programs entail in the compilation process. However, due to the inherent
constraints embedded within the R1CS form, this paper ultimately elects to use
data flow as a starting point for research.

Paper Organization The paper is organized as follows: In the next section, a
brief preliminary review of zero-knowledge proof and related tools. Section 3
provides the process of the proposed algorithm in this paper. The technical
exposition in Section 4 explains in detail the logic of the critical steps and
their formal description. Section 5 presents the specific categories of
benchmarks and their corresponding experimental results. Lastly, Section 6
concludes the present study.

##  II Preliminaries

This section introduces the basic concepts and principles of zero-knowledge
proof and discusses the role of R1CS and Circom in zero-knowledge-proof
systems. It also explores the limitations of existing normalization techniques
for R1CS and presents the proposed data-flow-based normalization generation
algorithm, which is motivated by these limitations.

###  II-A ZK-SNARKS

Zk-SNARK, which stands for Zero-Knowledge Succinct Non-interactive Argument of
Knowledge, is a type of zero-knowledge proof introduced in a 2014 paper [12].
The objective of zk-SNARK is to enable one party to prove to another that they
possess specific knowledge without revealing the knowledge itself, and to do
so concisely and efficiently. The working principle of Zk-SNARK can be
simplified into several steps. First, users convert the information they wish
to verify into a mathematical problem called a computation. This computation
can be implemented using any high-level programming language, such as OCaml,
C++, Rust, or hardware description languages like Circom.

Figure 1: The pipeline diagram for zk-SNARKS.

Next, the computation result is usually transformed into an arithmetic
circuit, a computational model for computing polynomials. An arithmetic
circuit consists of inputs and multiple gates, performing an essential
arithmetic operation such as addition or multiplication. The entire arithmetic
circuit can be used to generate a specific format (R1CS), which takes the form
of a constraint-based formula system (rank-1 constraint system) expressing the
constraints of the arithmetic circuit. A verifiable arithmetic circuit is one
of the inputs of the zk-SNARKs algorithm.

The Quadratic Arithmetic Program (QAP) is a variant of linear pcp and plays a
crucial role in the verifiable arithmetic circuit’s conversion into the QAP
format. QAP is a formula system that employs polynomials to represent the
behavior of the arithmetic circuit. By utilizing QAP, the security of the zk-
SNARKs algorithm is enhanced, along with the improvement of its implementation
efficiency [13].

Finally, there is the zk-SNARK stage, where the verifiable arithmetic circuit
and QAP are used to generate a proof. Zk-SNARKs are powerful privacy
protection protocols that can be utilized in digital payments, blockchain
technology, and other fields. They can verify the authenticity of information
while protecting the user’s privacy. Despite its relatively complex working
principle, zk-SNARK technology has found widespread application, bringing
higher security and privacy protection to the digital world.

###  II-B Circuit Language

ZKP technology can address several fundamental issues in the modern digital
world, including identity user’s verification without compromising their
private information and safeguarding privacy data from unauthorized
exploitation. Within ZKP, arithmetic circuits play a vital role in describing
and computing complex operations. These circuits consist of a series of logic
gates that can perform various arithmetic operations, such as addition,
multiplication, and division. By combining these basic operations, the complex
arithmetic circuits can be constructed to execute diverse computational tasks.

The circuit referred to here is a theoretical computational model, not an
actual electronic circuit.

This is the formal arithmetic circuit definition in theoretical computer
science [14].

######  Definition II.1.

A finite field field $\mathbb{F}$ is a field that contains a finite number of
integer elements.

$\mathbb{F}=\\{0,\ldots,p-1\\}\text{ for some prime }p>2$

The operations $+,\times,=$ on $\mathbb{F}$ should $\pmod{p}$ after
calculation.

######  Definition II.2.

A circuit is a triple $(M,L,G)$, where

  1. $\cdot$

$M$ is a set of values,

  2. $\cdot$

$L$ is a set of gate labels, each of which is a function from $M^{i}$ to $M$
for some non-negative integer $i$ (where $i$ represents the number of inputs
to the gate), and

  3. $\cdot$

$G$ is a labelled directed acyclic graph with labels from $L$.

######  Definition II.3.

A arithmetic circuit is a map $C:\mathbb{F}^{n}\rightarrow\mathbb{F}$, where
$\mathbb{F}$ is a finite field.

  1. 1.

It is a directed acyclic graph (DAG) where internal nodes are labeled
$+,-,\text{or }\times$ and inputs are labeled $1,x_{1},\ldots,x_{n}$, the
edges are wires or connections.

  2. 2.

It defines an n-variable polynomial with an evaluation recipe.

  3. 3.

Where $|C|=\\#\text{ gates in C}$.

This chapter mainly focuses on the Circom language, employed at the core step
of zk-SNARK protocols for describing arithmetic circuits. Within the framework
of ZKP systems, several commonly used arithmetic circuit description languages
exist, including Arithmetica, libsnark DSL, and Circom. These languages are
typically employed for building and verifying ZKP systems and can be used to
describe various arithmetic circuits, such as linear constraint systems (LCS),
bilinear pairings, and quadratic circuits.

By representing an arithmetic circuit as a constraint system, Circom’s core
idea involves describing inputs, outputs, and computation processes as linear
equations and inequalities. This approach allows developers to define complex
computation processes and generate the corresponding R1CS constraint system.
Circom provides a set of high-level abstract concepts, enabling developers to
focus more on the algorithm without being overwhelmed by low-level
implementation details.

Due to its comprehensive range of features, Circom has been widely used in
cryptography, blockchain, and other security-critical applications. Its
efficient and streamlined development of complex computational structures
gives developers a chance to quickly implement various privacy-preserving
protocols and zero-knowledge-proof techniques. Circom provides a simple,
declarative way of defining constraint systems and describing real-world
puzzles, which is of great importance to the protection of user’s privacy and
data security.

###  II-C Rank-1 Constraint Systems (R1CS)

R1CS, a common Arithmetic Circuit format that underlies real-world systems [2]
and an important part of the zk-SNARKS algorithm groth16[15], represents
computations as a set of constraint conditions, namely linear equations and
inequalities. Each equation has its own set of coefficients, while each
variable represents an input or output value. These equations and inequalities
describe the limiting conditions of the computation, implying that satisfying
these conditions correctly calculates the corresponding output result for the
given input sequence. R1CS includes a formal definition of constraint-based
computation rules, which can be verified using a set of public parameters and
a private input sequence. For a more detailed understanding of the formal
definition of R1CS, refer to Vitalik’s blog [13].

######  Definition II.4.

R1CS is a format for ZKP ACs. An R1CS is a conjunction of constraints, each of
the form:

$(\vec{a}\cdot\vec{x})\times(\vec{b}\cdot\vec{x})=(\vec{b}\cdot\vec{x})$

where $\vec{a},\vec{b}$ and $\vec{c}$ are vectors of coefficients (elements of
the prime field consisting of the integers modulo some prime), and $vec{x}$ is
a vector of distinct ”pseudo-variables”. Each pseudo-variable is either a
variable, representing an element of the field, or the special symbol 1,
representing the field element 1. $\cdot$ represents taking the dot product of
two vectors except that all additions and multiplications are done $\pmod{p}$,
and $\times$ represents the product of 2 scalars $\pmod{p}$. Using a pseudo-
variabale of 1 allows a constant addend to be represented in a dot product.

In Circom, each triplet of vectors in R1CS represents a mathematical
constraint. These vectors consist of coefficients of variables found at
corresponding positions in the solution vector $\vec{s}$. The solution vector
$\vec{s}$ includes assignments for all variables present in the R1CS
equations.

For example, a satisfied R1CS is shown in Fig.2:

Figure 2: A satisfied R1CS.

This constitutes a first-order constraint corresponding to a circuit
multiplication gate. If we combine all constraints, we obtain a first-order
constraint system.

R1CS is widely used in practical applications as a powerful computational
model. It serves as an integral component of the groth16 algorithm, which is a
popular version of zk-SNARK algorithms. R1CS plays a crucial role in improving
developers’ understanding of computer science and cryptography. Additionally,
it offers crucial support for various privacy protection measures.54

In this paper, we propose the R1CS paradigm for constraint groups. It imposes
constraints on the form and ordering of variable constraints.

######  Definition II.5.

R1CS paradigm is an R1CS satisfies the following requirements:

  1. 1.

If a constraint in the R1CS paradigm contains multiplication between
variables, it cannot have any other operators.

  2. 2.

If a constraint in the R1CS paradigm does not contain multiplication between
variables, it cannot contain intermediate variables generated by other linear
constraints.

  3. 3.

The ordering of constraints and variables (defined in Definition II.4) in the
R1CS paradigm must be consistent with the ordering method ($\times$ is in
front of $+$) in this paper.

The specific adjustments required for converting a general constraint system
into an R1CS paradigm will be explained through examples with concrete
constraints.

Requirement 1 suggests that complex quadratic constraints in the R1CS should
be split into simpler forms.

For instance.

$\displaystyle a\times b+c+d=f\Longrightarrow\ a\times b=r,r+c+d=f$
$\displaystyle 5\times a\times b=c\Longrightarrow 5\times a=r,r\times b=c$

Requirement 2 indicates that linear constraints in the R1CS system must be
eliminated by removing intermediate variables defined by other linear
constraints. For example,

$\displaystyle a+b=c,c+d=e\Longrightarrow a+b+d=e$

The specific sorting methods in requirement three will be discussed in later
sections outlining the algorithm’s steps.

###  II-D Data Flow Graph

In a bipartite-directed graph, known as a data flow graph, there are two types
of nodes: links and actors. Actors are utilized to represent various
operations, while links serve as the means by which data is received by
actors. Additionally, arcs allow links to transmit values to actors. The
formal definition of this concept can be found in Dennis’ paper [16].

######  Definition II.6.

A data flow graph is a bipartite labeled graph where the two types of nodes
are called actors and links.

$G=\langle A\cup L,E\rangle$ (1)

where

$\displaystyle A={a_{1},a_{2},\ldots,a_{n}}$ is the set of actors
$\displaystyle L={l_{1},l_{2},\ldots,l_{m}}$ is the set of links
$\displaystyle E\subseteq(A\times L)\cup(L\times A)$ $\displaystyle\text{is
the set of edges}.$

A more detailed description can be found in [17].

###  II-E Weighted Pagerank Algorithm

In this paper, we adopt the weighted PageRank algorithm to compute the weight
of each node in the data flow graph [18].

Pagerank algorithm is a method used for computing the ranking of web pages in
search engine results. It was initially proposed by Larry Page and Sergey
Brin, co-founders of Google, in 1998 and has since become one of the most
essential algorithms in the field of search engines.[19]

The algorithm assesses online web pages to determine their weight values,
which it then utilizes to rank search results. PageRank is based on the notion
that the weight of a web page is influenced by both the quantity and quality
of the other web pages that link to it.

The main steps of the Pagerank algorithm are as follows:

  1. 1.

Building the graph structure: First, the web pages and links on the internet
must be converted into a graph structure. In this structure, each web page
corresponds to a node and each link to a directed edge that points to the
linked web page.

  2. 2.

Computing the initial scores of each page: In Pagerank, the initial score of
each page is set to 1. This means that initially, each node has an equal
score.

Figure 3: The initial state of Pagerank algorithm.

  3. 3.

Iteratively computing the scores of each page: Each node’s score is
iteratively calculated based on its incoming links and averaged onto its
outgoing links at each iteration.

  4. 4.

Considering the number and quality of links: In addition to the relationships
between nodes, Pagerank considers the number and quality of links pointing to
a web page. Links from high-quality websites may carry more value than those
from low-quality sites. Therefore, when computing scores, the algorithm
weights links according to their number and quality.

Figure 4: The interation of Pagerank algorithm.

  5. 5.

Iterating until convergence: When the score of a node stabilizes, the
algorithm stops iterating. This indicates that the final scores of all nodes
have been determined and can be used to rank search results.

The Weighted PageRank algorithm differs from the standard PageRank algorithm
in that it incorporates the weight of each link as a factor, resulting in a
more precise evaluation of a webpage’s importance. Considering the importance
of pages, the original PageRank formula is modified as

$PR(u)=(1-d)+d\sum_{v\in B(u)}PR(v)W_{(v,u)}^{in}W_{(v,u)}^{out}$ (2)

In this equation, $W_{(v,u)}^{in}$ and $W_{(v,u)}^{out}$ are the weight of
$link(v,u)$ calculated based on the number of inlinks and outlinks of page $u$
and the number of inlinks of all reference pages of page $v$.

In this paper, we aim to use this algorithm to obtain more accurate weight
values for each node in the data flow graph.

sectionOverview In this section, we will introduce the procedure of
normalization through the process of converting R1CS introduced in Vitalik’s
blog [13] in the algorithm:

Constraint Set:

$\displaystyle A=\begin{pmatrix}0&1&0&0&0&0\\\ 0&0&0&1&0&0\\\ 0&1&0&0&1&0\\\
5&0&0&0&0&1\end{pmatrix}B=\begin{pmatrix}0&1&0&0&0&0\\\ 0&1&0&0&0&0\\\
1&0&0&0&0&0\\\ 1&0&0&0&0&0\end{pmatrix}$ $\displaystyle
C=\begin{pmatrix}0&0&0&1&0&0\\\ 0&0&0&0&1&0\\\ 0&0&0&0&0&1\\\
0&0&1&0&0&0\end{pmatrix}$

Firstly, the arithmetic tree generation process involves creating an
arithmetic tree for each constraint within the input R1CS constraint group,
which is subsequently merged. The resulting arithmetic tree comprises common
subformulas stored in a _DAG_. The constructed data flow graph’s structure is
shown in figure 6, which illustrates how the arithmetic trees are combined to
form the data flow graph.

Subsequently, a tile selection algorithm is implemented based on the data flow
graph, which divides the graph into tiles. The division of the entire graph
into tiles is illustrated in figure 6, depicting the overall procedure of the
tile selection algorithm. The specifics of the tile selection process,
including the form and selection logic, will be elaborated upon in subsequent
chapters.

Figure 5: The procedure of constructing the data flow graph.

Figure 6: The procedure of tile selection.

Next, the data flow graph is abstracted further with the selection of tiles as
a reference. A new abstracted node in the data flow graph replaces linear
constraints represented by tiles. The abstracted node can be represented as an
affine mapping, which preserves the linear relationship between the variables,
enabling faster computation of the intermediate values during the proof
generation process. This abstraction procedure streamlines the proof
generation process and reduces the computational cost of generating the proof.

Figure 7: The procedure of constructing the data flow graph.

We calculate the weight of each node with coefficients of the constraint. Then
we calculate the weights of the selected individual tiles using the improved
Weighted PageRank algorithm. The convergence process of the PageRank values of
four nodes in the abstract graph is depicted in figure 8.

Figure 8: The convergence process of each node.

Finally, constraints in the paradigm of R1CS are generated separately for each
tile. And the constraints and variables are ranked by the node weights
computed in the previous steps.

Now we convert the input R1CS to its paradigm:

$\displaystyle A=\begin{pmatrix}0&1&0&0&0\\\ 0&1&0&0&0\\\
5&1&0&1&-1\end{pmatrix}B=\begin{pmatrix}0&0&1&0&0\\\ 0&1&0&0&0\\\
1&0&0&0&0\end{pmatrix}$ $\displaystyle C=\begin{pmatrix}0&0&0&1&0\\\
0&0&1&0&0\\\ 0&0&0&0&0\end{pmatrix}$

##  III Normalization Algorithm

In this section, we formally introduce various steps of the normalization
generation algorithm and several data structures defined within the algorithm.

###  III-A Construction of RNode Graph

Our study presents a novel data structure, RNode, that represents variable
relationships within an R1CS arithmetic circuit. This structure, _RNode_ ,
facilitates efficient problem-solving by tracking interconnections among
variables, which is formally defined as follows,

######  Definition III.1.

An RNode is a node of two types in the data flow graph constructed in this
normalization algorithm.

$\displaystyle RNode$ $\displaystyle=ConstNode\cup VarNode$ (3) $\displaystyle
ConstNode$ $\displaystyle=\\{ConstValue,Operation,Father,Child\\}$
$\displaystyle VarNode$ $\displaystyle=\\{Operation,Father,Child\\}$

where

$\displaystyle\forall c\;\text{is a ConstNode}\cap
c.Operation=Null,c.child=\emptyset$ $\displaystyle\forall c\;\text{is a
ConstNode}\cap c.Operation\in\\{Add,Mul\\},c.father=\emptyset$

RNodes can be categorized into two types based on the variables they
represent. The first category represents the original variables in the
solution vector of the R1CS and the intermediate variables produced during the
construction of the arithmetic circuit. The second category represents the
constants in the data flow graph. Each RNode includes both an operator and a
computed result, which store the calculation method between its two parent
nodes and represent the calculated result of the subtree rooted in itself.

The generation of the RNode Graph involves 3 main stages:

  1. 1.

Transform each constraint into an equation of $a*b=c$, as required by the R1CS
constraints.

  2. 2.

Convert each constraint in the original R1CS constraint into an equation.

  3. 3.

Organize the resulting equations containing common sub-expressions into a DAG-
structured expression tree.

The core logic of the RNode Graph generation algorithm is similar to the
procedure we previously mentioned for constructing the RNode Graph.

The RNode is a data structure used to store information about the variables in
an R1CS during the construction of the RNode Graph. Unlike typical nodes in an
expression tree, each RNode stores information about both the variable and
operator involved in a given operation, allowing each operator’s output to be
considered an intermediate variable and making it more closely aligned with
the properties of an R1CS constraint set.

In practice, the merging and splitting of constraints poses a challenge in
determining the equivalence of R1CS constraint sets. However, our RNode graph
generation algorithm has observed that this merging or splitting does not lead
to substantial differences. When constraints are merged, one variable is
subtracted from the original constraint set. For instance, the merging process
can be exemplified by the following example.

$\begin{pmatrix}1&1&0&0\\\ 0&1&1&0\end{pmatrix}\cdot\begin{pmatrix}1&0&0&0\\\
1&0&0&0\end{pmatrix}=\begin{pmatrix}0&0&1&0\\\ 0&0&0&1\end{pmatrix}$
$\downarrow$
$\begin{pmatrix}1&2&0&0\end{pmatrix}\cdot\begin{pmatrix}1&0&0&0\end{pmatrix}=\begin{pmatrix}0&0&0&1\end{pmatrix}$
$\downarrow$
$\begin{pmatrix}1&2&0\end{pmatrix}\cdot\begin{pmatrix}1&0&0\end{pmatrix}=\begin{pmatrix}0&0&1\end{pmatrix}$

However, the subtracted variable will be added back into the RNode Graph as an
intermediate node in the sum-product expression during the construction of the
RNode Graph. The reverse is also true.

Further abstraction of the RNode Graph is required to eliminate the difference
observed in the graphs generated by equivalent R1CS due to the different
variable ordering in the constraint set. The main difference is observed in
the order of addition in constructing the expression of continuous addition.
At this stage, we do not have sufficient information to determine the
sequential execution order. During the algorithm execution flow, two variables
are randomly selected and added together, resulting in a different structure
in the graph.

###  III-B Tile Selection

Here, we categorize tiles into three types and give the formal definition of
tile.

######  Definition III.2.

Tile is a tree-like subgraph of the RNode Graph, representing a constraint in
R1CS.

  1. 1.

Quadratic: Tiles with the form $x*y=z$, where $x$, $y$, and $z$ are variables.

  2. 2.

MulLinear: Tiles whose root is obtained by multiplying its two parents, with
at least one of the parents being a constant.

  3. 3.

AddLinear: Tiles whose root is obtained by adding its two parents.

Tile is essentially a set of linear equations generated by applying certain
constraints to variables in the R1CS. We can use tiles as building blocks to
construct a normalized R1CS that is both correct and scalable.

During tile selection, we divide the data flow graph from the previous step
into 3 types. While AddLinear and MulLinear tiles are generated by linear
constraints and are essentially linear tiles, their logical processing differs
significantly. Hence, we discuss them as two separate types.

  1. 1.

We temporarily put aside the constraint merging step until we obtain more
information about the tree in subsequent steps.

  2. 2.

If there is a need to generate merged formulas later, it can be achieved
simply by applying a fixed algorithm to the unmerged formulas.

  3. 3.

The implementation of the tile selection algorithm is relatively simple.

When the RNode graph is generated, we select a node with no successors and use
it as the root to partition a subgraph from the RNode Graph as a tile. We then
remove this tile from the RNode Graph and repeat this loop until the entire
graph is partitioned. Detailed code will posted in github.

The difference between data flow graphs generated by equivalent R1CS
constraint systems lies in the order of node additions when processing linear
tiles. However, the nodes added within a linear constraint, when considered as
a set, are equivalent. This implies that the different addition order only
affects the traversal order of the nodes. Therefore, if we consider the
selected linear tiles as the products of the selected nodes and their
respective coefficients, the chosen sets of linear tiles from equivalent R1CS
constraint systems are the same. In other words, there is no distinction
between the selected sets of tiles for equivalent R1CS constraint systems.
Revised paragraph: The difference between data flow graphs generated by
equivalent R1CS constraint systems lies in the order of node additions when
processing linear tiles. However, the nodes added within a linear constraint,
when considered as a set, are equivalent. This implies that the different
addition order only affects the traversal order of the nodes. Therefore, the
selected sets of linear tiles from equivalent R1CS constraint systems are the
same if we consider the selected linear tiles as the products of the selected
nodes and their respective coefficients. In other words, there is no
distinction between the selected sets of tiles for equivalent R1CS constraint
systems.

###  III-C Graph Abstraction

We further abstract the data flow graph based on the selected tiles to
eliminate the differences between various equivalent R1CS constraint systems.
Specifically, we abstract linear tiles using a new node. By doing so, we mask
the difference in addition order within linear tiles in the RNode Graph and
transform the relationship between external nodes and specific nodes within
the linear tile into a relationship between external nodes and the tile to
which the particular node belongs. In this abstracted data flow graph, the
types of edges are as follows:

  1. 1.

Non-linear tile abstract node to non-linear tile abstract node: The two
vertices already existed in the original RNode graph. This edge type remains
consistent with the original RNode graph.

  2. 2.

Non-linear tile abstract node to linear tile abstract node: This edge type
exists only if there are non-abstract nodes in the linear tile represented by
the abstract node.

  3. 3.

Linear tile abstract node to linear tile abstract node: This edge type exists
only if the two abstract nodes represent linear tiles that share common non-
abstract nodes.

###  III-D Tile Weight Calculation

In this step, we use the Weighted PageRank algorithm to calculate the scores
of each vertex in the data flow graph.

The previous steps eliminated the differences in the data flow graphs
generated by equivalent R1CS through the abstraction of linear tiles. In the
following step, constraints are generated on a tile-by-tile basis, and a
criterion for tile order is proposed to sort the generated constraints.

In the algorithm proposed in this paper, the Weighted PageRank algorithm is
used to calculate the weights of each node in the abstracted data flow graph,
which are then used as the basis for calculating the weights of the
corresponding constraints for each tile. Compared to the traditional PageRank
algorithm, this algorithm assigns weights to every edge in the graph and
adjusts the iterative formula for node weights. In the Weighted PageRank
algorithm, as mentioned in the former section, the formula for calculating
node scores is defined as follows:

$PR(u)=(1-d)+d\sum_{v\in B(u)}PR(v)W_{(v,u)}^{in}W_{(v,u)}^{out}$ (4)

In this algorithm, the primary purpose of using the Weighted PageRank
algorithm is to reduce the symmetry of the abstracted data flow graph. The
structure of the graph is significantly simplified in the previous step
through the simplification of linear tiles. However, some structures in the
graph still contain symmetric nodes. If general algorithms are used to
calculate the weights of nodes, these symmetric nodes may be assigned the same
weight, which can lead to problems in subsequent constraint generation and
sorting. To address this issue, the Weighted PageRank algorithm is employed to
calculate weights for different nodes. This helps increase the asymmetry of
the graph and minimize the occurrence of nodes with the same score.

In this algorithm, the iterative formula for scores in the Weighted PageRank
algorithm is further adjusted. The node weights retained in the original data
flow graph are set to 1, while for nodes abstracted from linear constraints,
their weights are calculated through a series of steps.

First, for a linear constraint:

$\sum^{n}_{i=1}a_{i}b_{i}=c$ (5)

Convert it to:

$\sum^{n}_{i=1}a_{i}b_{i}-c=0$ (6)

Finally, the variance of the normalized coefficients is utilized as the weight
of the abstract node representing the linear constraint.

$W=\frac{\sum_{i=1}^{n}\left(a_{i}-\frac{\sum_{i=1}^{n}a_{i}-1}{n+1}\right)^{2}+\left(-1-\frac{\sum_{i=1}^{n}a_{i}-1}{n+1}\right)^{2}}{\left(\frac{\sum_{i=1}^{n}a_{i}-1}{n+1}\right)^{2}}$
(7)

In this algorithm, the iterative formula for the node score in the Weighted
PageRank algorithm is given by:

$PR(u)=(1-d)+d\sum_{v\in B(u)}PR(v)W^{u}W^{v}$ (8)

The scores for each node are calculated to determine the ranking of each
constraint in the R1CS form.

###  III-E Constraint Generation

In this step, we generate constraints in descending order of weight, with
tiles as the unit, starting with quadratic constraints and followed by linear
constraints.

In this phase, a portion of the variable ordering was determined through the
quadratic tile. In the tripartite matrix group of the R1CS paradigm, the
three-row vectors corresponding to quadratic constraints each have only one
non-zero coefficient. Still, in matrices A and B, the two-row vectors
corresponding to the same constraint are equivalent in the constraint
representation, consistent with the commutativity of multiplication. This
results in two equivalent expressions for the same quadratic constraint.
Figure 9 illustrates an example of equivalent terms for a quadratic
constraint.

Figure 9: Two equivalent expressions of a quadratic constraint.

In the abstract data flow graph, all vertices representing variables appearing
in quadratic constraints are retained, making it possible to determine the
choice of non-zero coefficients in the corresponding row vectors of matrices A
and B in the tripartite matrix group of the R1CS paradigm based on the weights
of each variable. The variables with higher weights are assigned to matrix A
and given smaller indices in the variable mapping.

The sorting rules for the ordering of variables that appear in quadratic
constraints can be summarized as follows:

  1. 1.

Sort variables based on the highest weight value among all quadratic
constraints in which they appear, such that variables with higher weight
values have smaller indices in the variable mapping.

  2. 2.

For variables that appear in the same constraint and have the same highest
weight value sort them based on their scores in the Weighted PageRank
algorithm for the nodes in the data flow graph corresponding to the variables.
Variables with higher node scores will have smaller indices in the variable
mapping and be assigned to the corresponding row vector in matrix A of the
tripartite matrix group.

###  III-F Adjustment of Linear Constraints

At this stage, the partition and ordering of constraints within a constraint
group have been established, and the ordering of variables that appeared in
quadratic constraints has also been determined. However, it is necessary to
adjust the order of newly introduced variables in linear tiles in this step.
As the example below demonstrates, introducing multiple new variables within a
linear tile can result in disorderly sorting. This is because, in the previous
actions, the specific structures of linear tile constraints were abstracted to
eliminate differences in the RNode Graph.

$\displaystyle A=\begin{pmatrix}0&1&0&0&0&0\\\ 0&1&0&0&0&0\\\
5&1&0&1&-1&1\end{pmatrix}B=\begin{pmatrix}0&0&1&0&0&0\\\ 0&1&0&0&0&0\\\
1&0&0&0&0&0\end{pmatrix}$ $\displaystyle C=\begin{pmatrix}0&0&0&1&0&0\\\
0&0&1&0&0&0\\\ 0&0&0&0&0&0\end{pmatrix}$ $\displaystyle
A=\begin{pmatrix}0&1&0&0&0&0\\\ 0&1&0&0&0&0\\\
5&1&0&1&1&-1\end{pmatrix}B=\begin{pmatrix}0&0&1&0&0&0\\\ 0&1&0&0&0&0\\\
1&0&0&0&0&0\end{pmatrix}$ $\displaystyle C=\begin{pmatrix}0&0&0&1&0&0\\\
0&0&1&0&0&0\\\ 0&0&0&0&0&0\end{pmatrix}$

Therefore, we introduce a new method to order the newly introduced variables
in linear tiles. For each new variable introduced in a linear tile, its weight
is calculated as

$weight=\sum_{other\ linear\ tiles}\left|field*weight\ of\ linear\
tile\right|$ (9)

The new variables are sorted based on their weights. If the weights are the
same, the coefficients of the variable in its linear tile are considered for
comparison. The appearance of new variables in other linear tiles, to some
extent, reflects their importance in the entire constraint group.
Additionally, suppose certain new variables only appear in their constraints.
In that case, their weights will be zero, and their ordering will only affect
the constraints generated by their linear tile without changing the ordering
of other constraints. Therefore, they can be sorted in descending order based
on their coefficients alone.

##  IV Evaluation

In this section, we introduced the self-designed benchmark used in this paper.
We evaluated the effectiveness of the paradigm generation algorithm by
analyzing the test results and the intermediate outputs.

###  IV-A Benchmark Design

To evaluate the proposed algorithm in this paper, we implemented the entire
process of paradigm generation explained in the former section using Python to
verify its results.

Due to the lack of related research, this field has no comprehensive
benchmark. Therefore, we summarized some rules for generating equivalent R1CS
constraint groups based on the logic the mainstream Circom compiler used to
create R1CS and designed a more comprehensive benchmark based on these rules.
The benchmark includes the following main categories depending on the
reflected situation:

  1. 1.

Replacement of variable order in R1CS.

  2. 2.

Transformation of constraint order in R1CS.

  3. 3.

Introduction of multiple new variables in a single linear constraint in R1CS.

  4. 4.

Introduction of new variables in multiple linear constraints in R1CS, with
shared new variables.

  5. 5.

Merging and splitting of constraints in R1CS.

The different categories in the benchmark correspond to the other reasons for
generating equivalent R1CS. Each category contains 2-3 basic R1CS constraints.
To comprehensively test the robustness and correctness of the algorithm, 5-6
equivalent R1CS constraint groups are generated for each R1CS based on the
respective reasons. The equivalent constraint groups of each constraint group
are paired and inputted into the algorithm to verify whether the algorithm can
generate consistent and R1CS-compliant output results defined in definition
II.5, when processing different equivalent constraint groups.

The publicly available data set used in this study can be found at the
following GitHub repository:
https://github.com/Ash1sc/R1CS_normalization_benchmark. This repository
contains the raw data that was utilized for testing purposes. It is important
to note that the data set is licensed under the GNU General Public License
version 3.0 (GPL-3.0). This license allows for the data set and scripts to be
freely distributed, modified, and used, with the condition that any derived
works are also licensed under the GPL-3.0 and that the original copyright and
license information is retained. If you would like to learn more about the
details of the GPL-3.0, please visit
https://www.gnu.org/licenses/gpl-3.0.en.html.

###  IV-B Result Evaluation

Table I shows the result of the experiments.

TABLE I: Experimental Results of Equivalent R1CS Constraint Group Conversion Reasons for Generating Equivalent R1CS Constraints |  Number of Groups |  Successfully Generated Groups | Pass Rate
---|---|---|---
Replacement of variable order in R1CS. |  55 |  55 | 100%
Reordering of constraint sequences in R1CS. |  21 |  21 | 100%
Introduction of multiple new variables in a single linear constraint in R1CS. |  15 |  15 | 100%
Introduction of multiple new variables with shared usage in multiple linear constraints in R1CS. |  15 |  15 | 100%
Merging and splitting of constraints in R1CS. |  6 |  6 | 100%

Observation shows that the generated paradigms meet the requirements of the
R1CS paradigm mentioned above and have equal semantics as the original R1CS
constraint groups.

Through analysis of the intermediate outputs at each stage of the conversion
process, it was found that for equivalent R1CS constraint groups generated by
reordering constraints, the only difference in the resulting data flow graphs
lies in the order in which nodes representing intermediate variables are
created. This is due to different processing orders of each constraint during
traversal, leading to other orders of introducing intermediate variables in
each constraint.

For equivalent R1CS constraint groups generated by variable replacement, the
differences lie in the order in which RNodes representing initial variables in
the R1CS are developed and in the order of addition in the summation chain
structure caused by differences in variable order in linear constraints.
However, these changes do not affect the selected tile set, and the same data
flow graph is obtained after abstraction.

In the final step of the algorithm, we proposed a novel variable ordering
method to solve the sorting confusion issue when multiple variables are
introduced to a linear constraint.Experimental results demonstrate that our
algorithm is capable of correctly identifying these variables and produces a
variable mapping sequence that conforms to the definition.

After the abstraction of the data flow graph, the splitting and merging of
linear constraints can lead to changes in the order of addition in the
summation chain. However, this issue is resolved. When equivalent R1CS
constraint groups are created through the merging and splitting of
constraints, the only discrepancy in the resulting data flow graphs is with
regards to the vertices representing intermediate variables. These vertices
can either represent existing variables or intermediate variables introduced
during the data flow graph’s creation. Nevertheless, the structure of the data
flow graph remains unaffected.

##  V Conclusion

In this paper, we propose an algorithm based on data flow analysis to
construct a paradigm for R1CS. The correctness and equivalence of R1CS have
long been challenging to study due to the diversity and flexibility of
constraint construction methods. Our algorithm aims to eliminate the
differences between equivalent R1CS constraint systems through a series of
abstraction processes. We introduce ordering methods for internal variables
and constraints in R1CS, which provide reference information for the final
paradigm generation. Through experimental results, we demonstrate that our
algorithm can identify equivalent R1CS resulting from constraint merging,
intermediate variable selection, and variable and constraint reordering, and
transform them into a unique paradigm. Therefore, R1CS plays an indispensable
part in zk-SNARKS. Further work includes establishing rules for merging
constraints and a more comprehensive benchmark. This requires us to conduct
more in-depth research and exploration of the generation rules of R1CS to
improve the algorithm.

## References

  * [1] Shafi Goldwasser, Silvio Micali and Chales Rackoff  “The knowledge complexity of interactive proof-systems”  In _Providing Sound Foundations for Cryptography: On the Work of Shafi Goldwasser and Silvio Micali_ , 2019, pp. 203–225
  * [2] “ZCash Company”, https://z.cash/, 2014
  * [3] Eli Ben Sasson et al.  “Zerocash: Decentralized anonymous payments from bitcoin”  In _2014 IEEE symposium on security and privacy_ , 2014, pp. 459–474  IEEE
  * [4] Guillaume Drevon  “J-R1CS – a JSON Lines format for R1CS”  In _Sikoba Research ∗_ , 2019
  * [5] Elvira Albert et al.  “Distilling constraints in zero-knowledge protocols”  In _Computer Aided Verification: 34th International Conference, CAV 2022, Haifa, Israel, August 7–10, 2022, Proceedings, Part I_ , 2022, pp. 430–443  Springer
  * [6] Eli Ben-Sasson et al.  “Aurora: Transparent succinct arguments for R1CS”  In _Advances in Cryptology–EUROCRYPT 2019: 38th Annual International Conference on the Theory and Applications of Cryptographic Techniques, Darmstadt, Germany, May 19–23, 2019, Proceedings, Part I 38_ , 2019, pp. 103–128  Springer
  * [7] Jonathan Lee, Srinath Setty, Justin Thaler and Riad Wahby  “Linear-time and post-quantum zero-knowledge SNARKs for R1CS”  In _Cryptology ePrint Archive_ , 2021
  * [8] Alexander Golovnev et al.  “Brakedown: Linear-time and post-quantum SNARKs for R1CS”  In _Cryptology ePrint Archive_ , 2021
  * [9] Yuan Muzi et al.  “Semantic comparison method and device between a kind of source code and binary code”, 2019
  * [10] Gao Li, Li Zhongqi, Yang Dongsheng and Yinzhong Liu  “Compiling method from intermediate language (IL) program to C language program of instruction list”  CN, 2013
  * [11] Zhao Yongsheng, Chen Zhiyong, Cui Rongtao and Wen Zhili  “A kind of assembly language is to the code conversion method of higher level language and device”  CN, 2015
  * [12] Eli Ben-Sasson, Alessandro Chiesa, Eran Tromer and Madars Virza  “Succinct $\\{$Non-Interactive$\\}$ zero knowledge for a von neumann architecture”  In _23rd USENIX Security Symposium (USENIX Security 14)_ , 2014, pp. 781–796
  * [13] Vitalik Buterin  “Quadratic arithmetic programs: from zero to hero”  In _URl: https://medium. com/@ VitalikButerin/quadratic-arithmetic-programs-from-zero-to-hero-f6d558cea649_ , 2016
  * [14] Heribert Vollmer  “Introduction to circuit complexity: a uniform approach”  Springer Science & Business Media, 1999
  * [15] Jens Groth  “On the size of pairing-based non-interactive arguments”  In _Advances in Cryptology–EUROCRYPT 2016: 35th Annual International Conference on the Theory and Applications of Cryptographic Techniques, Vienna, Austria, May 8-12, 2016, Proceedings, Part II 35_ , 2016, pp. 305–326  Springer
  * [16] Jack B Dennis  “First version of a data flow procedure language”  In _Programming Symposium: Proceedings, Colloque sur la Programmation Paris, April 9–11, 1974_ , 2005, pp. 362–376  Springer
  * [17] Philip C Treleaven, David R Brownbridge and Richard P Hopkins  “Data-driven and demand-driven computer architecture”  In _ACM Computing Surveys (CSUR)_ 14.1 ACM New York, NY, USA, 1982, pp. 93–143
  * [18] Wenpu Xing and Ali Ghorbani  “Weighted pagerank algorithm”  In _Proceedings. Second Annual Conference on Communication Networks and Services Research, 2004._ , 2004, pp. 305–314  IEEE
  * [19] Lawrence Page, Sergey Brin, Rajeev Motwani and Terry Winograd  “The pagerank citation ranking: Bring order to the web”, 1998


# Electronic Wigner-Molecule Polymeric Chains in Elongated Silicon Quantum
Dots and Finite-Length Quantum Wires

Arnon Goldberg<EMAIL_ADDRESS>Constantine Yannouleas
<EMAIL_ADDRESS>Uzi Landman
<EMAIL_ADDRESS>School of Physics, Georgia Institute of
Technology, Atlanta, Georgia 30332-0430

(03 February 2024)

###### Abstract

The spectral properties of electrons confined in a wire-like quasi-one-
dimensional (1D) elongated quantum dot (EQD) coupler between silicon qubits,
are investigated with a newly developed valley-augmented unrestricted Hartree-
Fock (va-UHF) method, generalized to include the valley degree of freedom
treated as an isospin, allowing calculations for a large number of electrons.
The lower energy symmetry-broken solutions of the self-consistent generalized
Pople-Nesbet equations exhibit, for a confinement that has been modeled after
an experimentally fabricated one in silicon, formation of Wigner-molecular
polymeric (longitudinal) chains, initiating through charge accumulation at the
edges of the finite-length quasi-1D wire. An increasing number of parallel
zig-zag chains form as the number of electrons loaded into the confinement is
increased, with the formation of newly added chains determined by the strength
of the transverse harmonic confinement. The broken-symmetry va-UHF solutions,
subsequently augmented by the quantum-mechanically required parity-
restoration, go beyond the va-UHF single-determinant solution, predicting
formation of entangled Wigner-molecular chains whose charge distributions
obliterate the zig-zag organization of the broken-symmetry solutions. The
symmetry-restored va-UHF methodology enables systematic investigations of
multi-electron complex nano-scale confined structures that could be targeted
for future imaging microscopy experiments in silicon and other materials
(e.g., 1D domain walls in TMD materials), and quantum information utilization.

##  I Introduction

Quantum-dot (QD) qubits are fundamental elements for semiconductor-based
solid-state quantum computing architectures [1, 2, 3]. One of the central
challenging issues in constructing scalable quantum processors is that of
quantum chip large-scale integration, allowing transfer of information between
computing qubits while preserving information during transfer.

Currently, attention focuses on patterned, gate-controlled elongated quantum
dots (EQDs) enabling coherent transfer of spins between relatively distant
quantum-dot qubits [4, 5, 6, 7], thus reducing the complexity and technical
difficulties that accompany short distance multi-dot couplers. To overcome the
challenges of designing qubit coupling strategies, recent developments [6, 7]
focus on silicon-based nanodevices [1, 2] guided by: (i) the long coherence
gained via the use of enriched 28Si substrates, (ii) the successful
demonstrations of Si-based high-speed operation and high fidelity spin quantum
qubits, and (iii) the vast investment made in the industry and the ensuing
infrastructure availability, as well as the scientific experience gained
already with Si-based technologies.

Here we aim at gaining fundamental insights about the many-body quantum nature
of the electronic states in such patterned, long-distance-coupling EQDs. Such
understanding is imperative for enabling theory-guided fabrication and
integration of these elements into solid-state Si-based quantum information
devices. Our main finding is that the extended nature of the (wire-like) EQD
results in conditions where the inter-electron repulsion energy dominates over
the electron quantal kinetic energy. These conditions underlie the emergent
formation of pinned Wigner Molecules (WMs), exhibiting a general architecture
of parallel chains. Quantum WMs (manifesting sliding 111In the case of a
confinement with circular symmetry, the WMs are referred to as rotating WMs as
well or pinned geometrical configurations) have been predicted theoretically
(see, e.g., Refs. [9, 10, 11, 12, 13, 14, 15, 16, 17]), and subsequently
confirmed experimentally in several materials systems (Ga[Al]As single QDs
[18, 19], GaAs double QDs [20, 21, 22], Si QDs [23, 24], carbon nanotubes
[25]), and most recently, in agreement with the latest theoretical predictions
[26, 27], in moiré superstructures in transition metal dichalcogenide (TMD)
materials [28].

Naturally, one proceeds by formulating the many-body Hamiltonian describing
the wire-like EQD, and by numerically solving the corresponding Schrödinger
equation. The most accurate method to this effect is exact diagonalization
(EXD) through the use of full configuration interaction (FCI) [29, 30, 14, 16,
22, 24], which is limited to systems containing up to about 10 confined
electrons, whereas the systems of relevance for qubit couplers are expected to
consist of a much larger number of electrons. Accordingly, we employ the
unrestricted Hartree-Fock (UHF) [31, 9, 11, 16], which refers to a family of
self-consistent field calculation techniques.

The UHF methodology uncovers formation of crystalline-like space-symmetry-
broken charge densities (CDs), which are finite-size analogs of the Wigner-
crystal chains investigated previously theoretically in the classical limit
[32] and in the context of Heisenberg chain modeling [33], as well as
investigated experimentally for very long GaAs/AlGaAs wires [34, 35].

A proper quantum mechanical description requires the restoration [11, 16, 36]
of the UHF broken symmetries with respect to all the symmetry operations
commuting with the system Hamiltonian. Such symmetry restoration is
implemented here with respect to the $y$-parity symmetry about the long
$x$-axis of the EQD. We discovered that the parity-symmetry-restored solutions
yield a substantially improved description of the CDs compared to those given
by the broken-symmetry UHF solutions. Moreover, this improvement, involving an
extension beyond the variational UHF single-determinant ansatz, proves to be
most efficient, entailing a rather minimal computational cost for systems
involving a relatively large number (dozens) of electrons.

Furthermore, the present implementation of the UHF method requires
consideration of the valley degree of freedom which is present in silicon
nanostructures, treated here as a pseudospin (along the spin). This requires
modifications (employed in this paper; see Sec. III) of the UHF Pople-Nesbet
equations [9, 11, 16] formulated originally in quantum chemistry for the case
of natural molecules [31]. We will interchangeably refer to our calculations
as UHF or valley-augmented UHF (va-UHF).

Figure 1: An illustration of the longboat-Viking-type total confining
potential $U(x,y)$ used. The square-type lateral potential, $U_{L}(x)$, along
the $x$-axis is given by Eq. (1). The transverse harmonic potential,
$U_{T}(y)$, along the $y$-direction corresponds to $\hbar\omega_{y}=10$ meV
and to an effective electron mass $m^{*}=0.19m_{e}$ (appropriate for silicon).
Figure 2: An illustration of the single-particle spectrum of electrons in the
potential, $U(x,y)$, used to model the elongated QD investigated in this
paper. For the parameters, see the text. The labels for the single-particle
space orbitals are as follows: $n_{x}\rightarrow$ number of nodes in the
$x$-direction. $n_{y}\rightarrow$ number of nodes in the $y$-direction.

##  II Many-body Hamiltonian

We consider a square-like, slightly asymmetric, confining potential along the
lateral long axis (the $x$-axis) of the EQD specified by:

$\displaystyle
U_{L}(x)=\xi|x/x_{0}|^{9}\Theta(-x)+\xi(x/x_{0})^{10}\Theta(x),$ (1)

where $\Theta(x)$ is the Heaviside step function, $x_{0}=50$ nm, and $\xi=0.1$
meV. This confinement corresponds to the EQD (referred to also as jellybean
QD) fabricated (including the above-noted slight asymmetry) and investigated
in Ref. [6]. The total confinement is given by

$\displaystyle U(x,y)=U_{L}(x)+U_{T}(y),$ (2)

where $U_{T}(y)=m^{*}\omega_{y}^{2}y^{2}/2$ is a harmonic potential along the
transverse $y$-direction; $m^{*}=0.19m_{e}$ is the effective electron mass for
Si and $\hbar\omega_{y}$ the harmonic quantum in the transverse $y$-direction.

An illustration of $U(x,y)$ is presented in Fig. 1. Furthermore, an
illustration of the associated (single-particle) spectrum (with
$\hbar\omega_{y}=2.5$ meV and $m^{*}=0.19m_{e}$) is plotted in Fig. 2.

The relevant many-body Hamiltonian is given by:

$\displaystyle H_{\rm MB}=\sum_{i=1}^{N}\left(\frac{{\bf
p}_{i}^{2}}{2m^{*}}+U({\bf
r}_{i})\right)+\frac{e^{2}}{\kappa}\sum_{i=1}^{N-1}\sum_{j>i}^{N}\frac{1}{|{\bf
r}_{i}-{\bf r}_{j}|},$ (3)

where $\kappa$ is the dielectric constant.

Figure 3: UHF ground-state charge densities (CDs) for $N=4-18$ electrons in a
silicon elongated quantum dot ($\sim 110$ nm long), with a transverse
confinemnet of $\hbar\omega_{y}=2.5$ meV. For each $N$ the CD is shown for
strongly repelling electrons ($\kappa=11$) at the top, and the CD for
electrons in the non-interacting limit is shown at the bottom. The non-
interacting CDs (at maximum occupation per orbital, large $\kappa=300$, and
$\hbar\omega_{y}=2.5$ meV) exhibit for all $N$ single-row structures of a
delocalized-particle nature, and obey a shell-filling Aufbau rule (following
spin and isospin exclusion rules; see text). In contrast, the strongly
interacting electrons CDs exhibit a transition from linear Wigner-molecular
structures for $N\leq 5$ to a single zig-zag chain at $N^{*}=5$ and to a
double zig-zag chain at $N^{*}=11$, that initiate via charge accumulation
occurring at the ends of the EQD. The color bars (bottom right) indicate the
charge density scale (in units of 1/nm2). For details, see the text.  Figure
4: UHF ground-state CDs for $N=1-20$ electrons for a silicon EQD with a
transverse confinement of $\hbar\omega_{y}=10$ meV. $\kappa=11$ appropriate
for silicon (strong interaction). The CDs exhibit a transition from one-row
Wigner-molecular structures for $N\leq 11$ to a zig-zag chain for $N\geq
N^{*}=12$ that initiates via charge accumulation occurring at the two ends of
the EQD (see $N=12-14$), with the zig-zag pattern fully developed for
$N=15-20$. For larger number of electrons a three-row zig-zag chain structure
emerges (not shown). The color bars (on the right) indicate the charge-density
scale (in units of 1/nm2).

##  III THE VALLEY-AUGMENTED UNRESTRICTED HARTREE-FOCK FOR NANOSYSTEMS

Due to strain in Si/SiGe quantum wells and (interfacial) quantum dots (in
particular, in heterostructure semiconductors) and higher subband quantization
energy in metal-oxide-semiconductor devices, the energies of the (four) in-
plane Si valleys are raised, resulting in two remaining degenerate valleys. In
this paper, as is the case for the general practice for the Si QDs [37], we
consider that the valley degree of freedom (VDOF) consists only of the low-
energy two-fold band.

As elaborated in Ref. [24], to characterize and classify the VDOF of the two
remaining valleys, we consider in this paper an isospin designation that is
constructed in analogy [including the SU(2) algebra generated by the
$i$-multiplied Pauli matrices] with that of the regular spin of the electrons
(replacing $\hat{\bf S}$ with $\hat{\bf V}$ when referring to the VDOF).
Obviously, in the absence of two-body interactions, the occupation of the
single-particle energy states 222These confinement-induced single-particle
energy states are also referred to as “orbitals”, see, e.g., the expressions
“atomic orbitals”, “space orbitals”, and “spin-orbitals” in chemistry [31] in
the QDs would depend on an interplay between the confinement (including an
applied magnetic field) and the spin and valley effects. For the case of two-
valley degeneracy (considered here) this interplay results in “quadrupling” of
the spectrum (four degenerate states for each confinement state, each
corresponding to a different spin and valley).

Here, we introduce the va-UHF for solving the many-body problem specified by
the Hamiltonian $H_{\rm MB}$ defined in Eq. (3). Accordingly, one introduces a
basis set of spatial single-particle orbitals $\varphi_{\mu}(x,y)$,
$\mu=1,\ldots,K$, which are given by the $K$ lowest-energy solutions of the
auxiliary single-particle Hamiltonian,

$\displaystyle H^{\rm core}=\sum_{i=1}^{N}\left(\frac{{\bf
p}_{i}^{2}}{2m^{*}}+U({\bf r}_{i})\right).$ (4)

For clarity and convenience, these solutions are sorted in ascending energy.

Subsequently, for the four families of va-UHF spin-isospin-orbitals 333This is
an apparent generalization of the term spin-orbital used in chemistry and
molecular physics, $\psi_{i}^{\alpha\zeta}\alpha\zeta$,
$\psi_{i}^{\beta\zeta}\beta\zeta$, $\psi_{i}^{\alpha\eta}\alpha\eta$,
$\psi_{i}^{\beta\eta}\beta\eta$, where $\alpha(\beta)$ denote up (down) spins,
$\zeta(\eta)$ denote up (down) isospins [i.e., electrons in the first (second)
valley], one considers the expansions:

$\displaystyle\begin{split}&\psi_{i}^{\alpha\zeta}=\sum_{\mu=1}^{K}C_{\mu
i}^{\alpha\zeta}\varphi_{\mu},\\\
&\psi_{i}^{\beta\zeta}=\sum_{\mu=1}^{K}C_{\mu i}^{\beta\zeta}\varphi_{\mu},\\\
&\psi_{i}^{\alpha\eta}=\sum_{\mu=1}^{K}C_{\mu i}^{\alpha\eta}\varphi_{\mu},\\\
&\psi_{i}^{\beta\eta}=\sum_{\mu=1}^{K}C_{\mu i}^{\beta\eta}\varphi_{\mu},\\\
&i=1,2,\ldots,K.\end{split}\ $ (5)

Then, following similar steps as in Ch. 3.8.2 of Ref. [31], we have derived a
generalization of the Pople-Nesbet equations as follows:

$\displaystyle\begin{split}&\sum_{\nu}F_{\mu\nu}^{\alpha\zeta}C_{\nu
j}^{\alpha\zeta}=\varepsilon_{j}^{\alpha\zeta}\sum_{\nu}S_{\mu\nu}C_{\nu
j}^{\alpha\zeta},\\\ &\sum_{\nu}F_{\mu\nu}^{\beta\zeta}C_{\nu
j}^{\beta\zeta}=\varepsilon_{j}^{\beta\zeta}\sum_{\nu}S_{\mu\nu}C_{\nu
j}^{\beta\zeta},\\\ &\sum_{\nu}F_{\mu\nu}^{\alpha\eta}C_{\nu
j}^{\alpha\eta}=\varepsilon_{j}^{\alpha\eta}\sum_{\nu}S_{\mu\nu}C_{\nu
j}^{\alpha\eta},\\\ &\sum_{\nu}F_{\mu\nu}^{\beta\eta}C_{\nu
j}^{\beta\eta}=\varepsilon_{j}^{\beta\eta}\sum_{\nu}S_{\mu\nu}C_{\nu
j}^{\beta\eta},\\\ &j=1,2,\ldots,K.\end{split}$ (6)

with the Fock-operator matrices being given by,

$\displaystyle\begin{split}F_{\mu\nu}^{\alpha\zeta}=H_{\mu\nu}^{\rm
core}+\sum_{\lambda}\sum_{\sigma}P_{\lambda\sigma}^{T}[(\mu\nu|\sigma\lambda)-P_{\lambda\sigma}^{\alpha\zeta}(\mu\lambda|\sigma\nu),\\\
F_{\mu\nu}^{\beta\zeta}=H_{\mu\nu}^{\rm
core}+\sum_{\lambda}\sum_{\sigma}P_{\lambda\sigma}^{T}[(\mu\nu|\sigma\lambda)-P_{\lambda\sigma}^{\beta\zeta}(\mu\lambda|\sigma\nu),\\\
F_{\mu\nu}^{\alpha\eta}=H_{\mu\nu}^{\rm
core}+\sum_{\lambda}\sum_{\sigma}P_{\lambda\sigma}^{T}[(\mu\nu|\sigma\lambda)-P_{\lambda\sigma}^{\alpha\eta}(\mu\lambda|\sigma\nu),\\\
F_{\mu\nu}^{\beta\eta}=H_{\mu\nu}^{\rm
core}+\sum_{\lambda}\sum_{\sigma}P_{\lambda\sigma}^{T}[(\mu\nu|\sigma\lambda)-P_{\lambda\sigma}^{\beta\eta}(\mu\lambda|\sigma\nu),\end{split}$
(7)

and the $\varepsilon_{j}^{\alpha\zeta}$, $\varepsilon_{j}^{\beta\zeta}$,
$\varepsilon_{j}^{\alpha\eta}$, $\varepsilon_{j}^{\beta\eta}$ being the
energies for the spin-isospin va-UHF orbitals. Because the $\varphi_{\mu}$’s
are eigenfunctions of $H^{\rm core}$, one has $S_{\mu\nu}=\delta_{\mu\nu}$ for
their overlaps. The two-body Coulomb matrix elements are given by

$\displaystyle\begin{split}(\mu\nu|\sigma\lambda)=\frac{e^{2}}{\kappa}\int
d{\bf r}_{1}d{\bf r}_{2}\varphi^{*}_{\mu}({\bf r}_{1})\varphi_{\nu}({\bf
r}_{1})r^{-1}_{12}\varphi^{*}_{\sigma}({\bf r}_{2})\varphi_{\lambda}({\bf
r}_{2}),\end{split}$ (8)

with $r^{-1}_{12}=1/|{\bf r}_{1}-{\bf r}_{2}|$.

In Eq. (7), the partial density matrices are given by

$\displaystyle\begin{split}P_{\mu\nu}^{\alpha\zeta}=\sum_{a}^{N^{\alpha\zeta}}C_{\mu
a}^{\alpha\zeta}(C_{\nu
a}^{\alpha\zeta})^{*},\;P_{\mu\nu}^{\beta\zeta}=\sum_{a}^{N^{\beta\zeta}}C_{\mu
a}^{\beta\zeta}(C_{\nu a}^{\beta\zeta})^{*},\\\
P_{\mu\nu}^{\alpha\eta}=\sum_{a}^{N^{\alpha\eta}}C_{\mu a}^{\alpha\eta}(C_{\nu
a}^{\alpha\eta})^{*},\;P_{\mu\nu}^{\beta\eta}=\sum_{a}^{N^{\beta\eta}}C_{\mu
a}^{\beta\eta}(C_{\nu a}^{\beta\eta})^{*},\end{split}$ (9)

the total-density matrix is defined as

$\displaystyle\begin{split}{\bf P}^{T}={\bf P}^{\alpha\zeta}+{\bf
P}^{\beta\zeta}+{\bf P}^{\alpha\eta}+{\bf P}^{\beta\eta},\end{split}$ (10)

and

$\displaystyle\begin{split}N^{\alpha\zeta}+N^{\beta\zeta}+N^{\alpha\eta}+N^{\beta\eta}=N.\end{split}$
(11)

The energy eigenvalues ($\varepsilon_{j}^{\alpha\zeta}$,
$\varepsilon_{j}^{\beta\zeta}$, $\varepsilon_{j}^{\alpha\eta}$,
$\varepsilon_{j}^{\beta\eta}$) and expansion coefficients ($C_{\nu
j}^{\alpha\zeta}$, $C_{\nu j}^{\beta\zeta}$, $C_{\nu j}^{\alpha\eta}$, $C_{\nu
j}^{\beta\eta}$) are obtained via self-consistent solutions of Eqs. (6).

We note that, unlike the valleytronic FCI [24], the va-UHF does not conserve
the total spin $\hat{\bf S}^{2}$ and total valley $\hat{\bf V}^{2}$ quantum
numbers; it only conserves their projections $S_{z}$ and $V_{z}$.

##  IV va-UHF charge densities

For $N=4-18$, fully spin and valley polarized ($S_{z}=V_{z}=N/2$), electrons
and $\hbar\omega_{y}=2.5$ meV, Fig. 3 contrasts the UHF CDs at $\kappa=11$
(silicon, top rows) to those associated with the non-interacting limit (NIL)
(at maximum occupation per orbital, $\kappa=300$, and $\hbar\omega_{y}=2.5$
meV, bottom rows). For all the strongly-interacting-electron cases, the UHF
densities exhibit explicit configurations of $N$ well defined humps reflecting
formation of pinned WMs, associated with the regime of strong inter-electron
correlations. In contrast, the NIL CDs conform to those expected from (Aufbau)
shell closures in the $U(x,y)$ confinement, with a fourfold period (due to
both the spin and valley) being clearly visible for $N=5-8$ [when the
$(1_{x},0_{y})$ orbital of the confining potential $U(x,y)$ is sequentially
occupied (see also Fig. 6); here $(n_{x},n_{y})$ denotes the number of nodes
in the $x$\- and $y$-directions].

In Fig. 8 of Appendix A, we demonstrate that the CDs shown in Fig. 3 (see also
Fig. 9), calculated for the spin and valley fully polarized
($S_{z}=V_{z}=N/2$) electrons, are essentially identical to the CDs for the
va-UHF states with minimal spin, $S_{z}$, and isospin, $V_{z}$, projections (0
for even $N$ or $\pm$1/2 for odd $N$). This finding concurs with the magneto-
spectroscopy measurements [6] where it has been found that the “Jellybean”
quantum dot studied in these experiments lacks visible spin structure for
similar electron numbers as those investigated in the present paper.

We recall that the strength of correlations is often expressed via the Wigner
parameter [9, 16] $R_{W}={\cal E}_{C}/\delta$, defined as the ratio of a
typical inter-electron repulsive Coulombic energy, ${\cal E}_{C}$, over a
typical energy gap, $\delta$, in the single-particle spectrum of $U(x,y)$;
naturally, the strong-interaction regime is expected for $R_{W}>1$,
accompanied by WM formation. For the square-like confinement here, we take
${\cal E}_{C}$ as the two-body matrix element of the Coulomb interaction for
two electrons occupying the nodeless lowest $(0_{x},0_{y})$ orbital, and for
$\delta$ we take the difference between the energies of the $(0_{x},0_{y})$
and $(1_{x},0_{y})$ orbitals. For the EQD investigated here, we have
$R_{W}=23.04$ for $\kappa=11$ and $\hbar\omega_{y}=2.5$ meV, $R_{W}=29.14$ for
$\kappa=11$ and $\hbar\omega_{y}=10$ meV (strong interactions, case of the Si
QD), and $R_{W}=0.84$ for $\kappa=300$ and $\hbar\omega_{y}=2.5$ meV (non-
interacting limit).

A prominent feature of the WM configurations is the successive formation with
increasing $N$ of complex, polymeric, multiple chain-like arrangements,
reminiscent of the classical Wigner-crystalline chains associated with the
equilibrium configurations of classical point charges [32] in infinite-length
wires. Naturally, the presence of the square-like edges here perturbs the
perfect-chain formations of the infinite-length wires by forcing an
accumulation of charges at the edges. The well-known zig-zag configuration,
however, is visible away from the edges of the EQDs for $N=7-11$ (single zig-
zag) and for $N=14-18$ (double zig-zag).

We note that a single row appears for $N=1-5$, whilst a second row starts
developing at $N=7$ with $N=6$ being a transitional stage. A third row starts
developing at $N=11$, whilst transitional cases towards a fourth row appear at
$N=17$ and $N=18$, with four electrons accumulating in a line along the
$y$-direction at the edges.

Fig. 4 displays the UHF CDs for fully spin and valley polarized $N=1-20$
electrons for a tighter transverse confinement with $\hbar\omega_{y}=10$ meV.
These CDs exhibit again organization with Wigner-chain-like features, the main
difference from the $\hbar\omega_{y}=2.5$ meV case [see Fig. 3] being a delay
in the appearance (as a function of $N$) of the transition regions between the
multiple chains. In particular, the single chain transitions to a double zig-
zag chain at the region $N=12-14$, with the double chain starting to form from
the edges 444This behavior contrasts with that in harmonic confinements along
the lateral $x$-direction, where the zig-zag chain starts forming at the
center of the linear chain, a fact well known from classical calculations in
the literature of trapped heavy ions; see, e.g., Refs. [46, 44, 45], a
behavior related to the accumulation of charge at the sharp square-like edges
of the confinement along the $x$-direction. Note that, due to this delay, the
third row does not appear in Fig. 4 in the range $N=1-20$; for
$\hbar\omega_{y}=10$ meV, a third row is expected to develop at larger values
of $N$.

Figure 5: Parity-restored CDs for a silicon EQD. (a) transverse confinement of
$\hbar\omega_{y}=2.5$ meV. $N=8,9,11,14$ fully spin and valley polarized
electrons. (b) transverse confinement of $\hbar\omega_{y}=10$ meV. $N=15-18$
fully spin and valley polarized electrons. In both cases, $\kappa=11$
appropriate for silicon (strong interaction). The symmetry-broken UHF CDs,
corresponding to the restored ones shown in (a) and (b), are given in Figs. 3
and 4, respectively. Obviously, only symmetry-broken UHF wave-functions
undergo symmetry restoration. CDs in units of 1/nm2.  Figure 6: UHF second
energy differences, $\Delta_{2}$. The curves associated with
$\hbar\omega_{y}=2.5$ meV and $\hbar\omega_{y}=10$ meV for $\kappa=11$ (strong
interaction), as well as with the non-interacting limit (at maximum occupation
per orbital, $\kappa=300$, and $\hbar\omega_{y}=2.5$ meV) are as noted. The
arrows denote the transition sizes ($N^{*}$’s) between the multiple chains
marked in Fig. 3 and Fig. 4.

##  V Restoration of parity along the transverse $y$-direction

The broken-symmetry UHF solutions violate a fundamental axiom of quantum
mechanics, namely, that the single-particle CDs must preserve the symmetries
of the many-body Hamiltonian. Here we take a first step in rectifying this UHF
inadequacy by restoring the parity symmetry along the transverse
$y$-direction, which is associated with a symmetric harmonic confinement
$U_{T}(-y)=U_{T}(y)$; note that the confinement along $x$ was taken to be
slightly anisotropic.

To restore the $y$-parity symmetry, we apply on the UHF Slater determinant
(SD), $\Psi_{\rm UHF}$, the projection operator [16, 36, 41]

$\displaystyle\hat{\Pi}=\frac{1}{2}(1+p\hat{\cal P}_{y}),$ (12)

with $p=\pm 1$. $\hat{\cal P}_{y}$ is the many-body parity operator, which
inverts about the $x$-axis the $y$-coordinates for all electrons, namely,
$\hat{\cal P}_{y}=\prod_{i=1}^{N}\hat{P}^{(i)}_{y}$.

The charge density associated with the projected wave function $\Psi_{\rm
PUHF}=\hat{\Pi}\Psi_{\rm UHF}$ is given by

$\displaystyle\rho_{\rm PUHF}({\bf r})={\cal N}\langle\Psi_{\rm
UHF}|\hat{\Pi}\sum_{i=1}^{N}\delta({\bf r}_{i}-{\bf r})\hat{\Pi}|\Psi_{\rm
UHF}\rangle,$ (13)

where the factor ${\cal N}$ imposes normalization to $N$. The expansion of the
symmetry-restored charge density in Eq. (13) is explicitly displayed in Eq.
(15) of Appendix C.

The projected total energy is given by [16, 36]

$\displaystyle E_{\rm PUHF}=\frac{\langle\Psi_{\rm UHF}|H\hat{\Pi}|\Psi_{\rm
UHF}\rangle}{\langle\Psi_{\rm UHF}|\hat{\Pi}|\Psi_{\rm UHF}\rangle}.$ (14)

Eqs. (13) and (14) are calculated using the property
$\hat{P}^{(i)}_{y}Y_{n}y)=\varpi_{n}Y_{n}(y)$, where $\varpi_{n}=(-1)^{n}$,
$n=0,1,2,\ldots$ and $Y_{n}(y)$ are the good-parity eigenfunctions of
$U_{T}(y)$.

Several CD results of this parity restoration are displayed in Fig. 5; it is
apparent that the zig-zag motif gets obliterated by the parity restoration.
Whether a zig-zag or a symmetry-restored CD will be the actual finding in a
quasi-1D system may depend in practice on several additional factors,
including the presence of impurities, the length of the system, as well as the
value of the effective mass. Indeed, the longer the system, the higher the
expectation for a broken symmetry CD. Interestingly, because of the much
larger mass compared to electrons, the CDs of trapped heavy ions are routinely
found experimentally to exhibit broken-symmetry structures, including zig-zag
chains. Nevertheless, in spite of the large mass of the ions, the
superposition (or entanglement) of both the “zig” and “zag” configurations is
quantum mechanically allowable, and recently it has been observed
experimentally [42].

Figure 7: Comparison of (a) symmetry-broken UHF, (b) $y$-parity-restored UHF,
and (c) FCI CDs for the associated ground states of $N=4$ electrons at zero
magnetic field in a single elliptic quantum dot, with an effective mass
$m^{*}=0.067m_{e}$ and a potential confinement specified by
$\hbar\omega_{x}=3.00$ meV and $\hbar\omega_{y}=5.64$ meV. The dielectric
constant used was $\kappa=1$. In (a) and (b), $S_{z}=0$. In (c) $S=S_{z}=0$.
Lenghts in units of nm. Charge densities in units of 1/nm2. The charge
densities are normalized to the total number of fermions $N=4$. The dashed
lines are a guide to the eye; in (b) and (c), they delineate two mirror
trapezoids.

Even though it is expected that (due to the much smaller electronic mass) the
actual CDs of electrons in EQDs and finite length wires would conform to those
produced through the symmetry-restoration correction, a promising experimental
opportunity emerged recently in the case of 1D domain walls in TMD bilayers
[43], where rapid progress in STM imaging could differentiate between
symmetry-broken and symmetry-preserving CDs.

##  VI Second energy differences

The second difference of the total energies, $\Delta_{2}=E(N+1)+E(N-1)-2E(N)$,
corresponding to the va-UHF electronic configurations shown in Figs. 3 and 4
are shown in Fig. 6. Also shown is the non-interacting case (at maximum
occupation per orbital, $\kappa=300$, and $\hbar\omega_{y}=2.5$ meV), which
exhibits shell closures at $N=4$, 8, 12, and 16 according to the Aufbau
principle, with the period of 4 resulting from both the spin and valley
degrees of freedom. This non-interacting behavior contrasts sharply with that
of the strong-interaction cases, where polymeric WM chains are formed. We note
that $\Delta_{2}$ is proportional to the capacitance of the EQD.

##  VII AN EXAMPLE OF A COMPARISON WITH full configuration interaction

In Fig. 7, we present a case of a comparison between UHF and parity-restored
UHF charge densities with those from the corresponding FCI [22, 24]
calculation. We choose the case of an elliptic external 2D potential that
binds $N=4$ conduction electrons for a one-band (no valley present)
semiconductor material.

The UHF CD in Fig. 7(a) displays a broken-symmetry zig-zag configuration.
However, restoration of the parity along the $y$-axis in Fig. 7(b) (which also
restores the $x$-parity in this example) produces a configuration of two
mirror trapezoids, which is formed by the spreading out (along the $x$-axis)
of the two middle sharp single-electron humps in panel Fig. 7(a). In this
context, note the reduction in height of the corresponding broad humps in Fig.
7(b). The parity-restored UHF CD in Fig. 7(b) is in good agreement with the
FCI CD in Fig. 7(c). Naturally, the FCI solution displays a certain degree of
additional relaxation effects in the CD humps.

##  VIII Discussion and Conclusions

We introduced a generalization (referred to as va-UHF) of the unrestricted
Hartree-Fock methodology [31, 9, 11, 16] that incorporates the valley degree
of freedom on an equal footing with the electronic spin. Furthermore, we
explored an augmenting symmetry (here, parity)-restoration step applied to the
symmetry-broken va-UHF (going beyond the single determinant description) that
amplifies and improves the quantum mechanical description of the considered
nanosystem [11, 16, 36]. This new va-UHF methodology is able to effectively
and economically investigate the physics of double-quantum-dot-based Si qubits
for a number of charge carriers much larger than the number accessible with
FCI approaches [29, 30, 14, 16, 22, 24]. The close agreement between the
results of charge densities obtained via the symmetry-restored (SR)-UHF
methodology and those evaluated by full configuration (exact diagonalization)
computations is illustrated in Fig. 7).

Figure 8: Ground-state va-UHF charge densities for $N=4-9$ fermions in the
same elongated-dot confinement as employed in Fig. 3 of the main text,
associated with states with minimum spin and isospin polarizations, as
follows: (A) $S_{z}=0$, $V_{z}=0$, $N=4$. (B) $S_{z}=1/2$, $V_{z}=1/2$, $N=5$.
(C) $S_{z}=0$, $V_{z}=0$, $N=6$. (D) $S_{z}=1/2$, $V_{z}=-1/2$, $N=7$. (E)
$S_{z}=0$, $V_{z}=0$, $N=8$. (F) $S_{z}=1/2$, $V_{z}=1/2$, $N=9$. The square-
like potential $U_{L}(x)$ along the $x$-direction is given by Eq. (1) in the
main text with $x_{0}=50$ nm and $\xi=0.1$ meV. The harmonic confinement
$U_{T}(x)$ along the $y$-direction has $\hbar\omega_{y}=2.5$ meV. Other
parameters used: effective mass $m^{*}=0.19m_{e}$ and dielectric constant
$\kappa=11$. Charge densities in units of 1/nm2.  Figure 9: The fully spin and
isospin polarized, $(S_{z}=V_{z}=N/2)$, va-UHF CDs for $N=4-9$ to be
contrasted with those in Fig. 8. Parameters are the same as in the
corresponding frames in Fig. 8. CDs in units of 1/nm2.

A pioneering va-UHF application was presented here for the case of an
elongated [7] (referred to also as a jellybean [6]) Si QD. EQDs have been
proposed as long-distance couplers between solid-state qubits, enabling thus
an effective solution to the problem of scalability [4, 5, 6, 7]. Our
calculations for a typical Si EQD revealed a strongly correlated regime as a
particular case of Wigner molecularization [9, 11, 16, 17], with the electrons
organized in polymeric multi-chain arrangements which are finite-size analogs
of the infinite Wigner chains considered earlier [32, 33, 34, 35]. In this
context, we note that our results agree with the experimental findings of Ref.
[7] that the charge is "well distributed" along the long lateral axis of the
EQD, whereas Ref. [6] concluded that the electrons bunch together forming
several separated QDs. Instead of separated QDs, our investigations suggest
that the fragmentation in the transverse direction associated with the
formation of multiple chains may explain the variety of conductance behaviors
reported in Ref. [6]. Conductance calculations between quantum dots coupled to
the elongated coupler element, based on the many-body restored-parity va-UHF
wave functions, will constitute a promising follow up step of the present
paper.

Finally, we note that the va-UHF findings could be further confirmed with a
newly developing experimental platform associated with 1D domain walls in
twisted TMD bilayers [43].

##  IX acknowledgments

This work has been supported by a grant from the Air Force Office of
Scientific Research (AFOSR) under Grant No. FA9550-21-1-0198. Calculations
were carried out at the GATECH Center for Computational Materials Science.

Figure 10: The energy difference (top), $\Delta E_{\pm}=E_{\rm PUHF}-E_{\rm
UHF}$ and expectation value of the parity operator (bottom), $\langle{\rm
UHF}|{\cal P}|{\rm UHF}\rangle$, for $N=7$ fully spin and valley polarized
electrons as a function of the dielectric conctant $\kappa$. The $\pm$
subscript corresponds to the two eigenvalues $p=\pm 1$ entering in the
definition of the parity projection operator; see Eq. (12) of the main text.
The parameters used are: $m^{*}=0.19m_{e}$, $\hbar\omega_{y}=2.5$ meV.

##  Appendix A SOME CASES OF CHARGE DENSITIES FOR va-UHF STATES WITH MINIMUM
SPIN AND ISOSPIN POLARIZATION

In Fig. 8, we demonstrate that the va-UHF CDs shown in Fig. 3, calculated for
the spin and valley fully polarized electrons, are essentially identical to
the CDs for the va-UHF states with minimal spin, $S_{z}$, and isospin,
$V_{z}$, projections; see the caption for the specific values. This finding
concurs with the magneto-spectroscopy measurements [6] where it has been found
that the “Jellybean” quantum dot studied in these experiments lacks visible
spin structure for similar electron numbers as those investigated in the
present paper. In order to facilitate the direct comparison, we display in
Fig. 9 the corresponding fully spin and isospin polarized ($S_{z}=V_{z}=N/2$)
va-UHF CDs.

##  Appendix B AN EXAMPLE OF CALCULATIONS OF PROJECTED ENERGIES

According to Eqs. (12)-(14) of the main text, there are always two projected
states, $\Psi_{\rm PUHF}$, with good parity eigenvalues $p=\pm 1$ for each
mixed-parity state, $\Psi_{\rm UHF}$. The top panel of Fig. 10 demonstrates
that one of the symmetry-restored states has always an energy lower (at most
equal to) than the broken symmetry UHF state. The bottom panel of Fig. 10
displays the expectation value of the parity operator, $\langle{\rm UHF}|{\cal
P}|{\rm UHF}\rangle$, which varies between -1 and +1. Note that, as a function
of a increasing dielectric constant $\kappa$, the negative-parity projected
ground state is replaced by the positive-parity one, following closely a
similar behavior by the expectation value.

##  Appendix C CLARIFICATION CONCERNING THE COMPACT EXPRESSION IN EQ. (13) OF
THE MAIN TEXT

Denoting the va-UHF Slater determinant as $\Psi(x,y)$ (we drop the subscript
UHF here), its mirror image about the $x$ axis is given by $\Psi(x,-y)$, and
the $y$-parity restored wave function is $\propto\Psi(x,y)+p\Psi(x,-y)$, with
$p=\pm 1$. Then, because the Slater determinant $\Psi(x,-y)$ is in general not
orthogonal to $\Psi(x,y)$, the expectation value of an operator ${\cal O}$ is
given by

$\displaystyle\frac{\langle\Psi(x,y)|{\cal
O}|\Psi(x,y)\rangle+p\langle\Psi(x,y)|{\cal
O}|\Psi(x,-y)\rangle+p\langle\Psi(x,-y)|{\cal
O}|\Psi(x,y)\rangle+\langle\Psi(x,-y)|{\cal
O}|\Psi(x,-y)\rangle}{\langle\Psi(x,y)|\Psi(x,y)\rangle+p\langle\Psi(x,y)|\Psi(x,-y)\rangle+p\langle\Psi(x,-y)|\Psi(x,y)\rangle+\langle\Psi(x,-y)|\Psi(x,-y)\rangle}.$
(15)

The operator associated with the charge density is a one-body operator,
$\sum_{i=1}^{N}\delta({\bf r}-{\bf r}_{i})$.

## References

  * Zwanenburg  _et al._ [2013a] F. A. Zwanenburg, A. S. Dzurak, A. Morello, M. Y. Simmons, L. C. L. Hollenberg, G. Klimeck, S. Rogge, S. N. Coppersmith, and M. A. Eriksson, Silicon quantum electronics, Rev. Mod. Phys. 85, 961–1019 (2013a).
  * Warren and Economou [2022] A. Warren and S. E. Economou, Silicon qubits move a step closer to achieving error correction, Nature 601, 320–322 (2022).
  * Burkard  _et al._ [2023] G. Burkard, T. D. Ladd, A. Pan, J. M. Nichol, and J. R. Petta, Semiconductor spin qubits, Rev. Mod. Phys. 95, 025003 (2023).
  * Martins  _et al._ [2017] F. Martins, F. K. Malinowski, P. D. Nissen, S. Fallahi, G. C. Gardner, M. J. Manfra, C. M. Marcus, and F. Kuemmeth, Negative Spin Exchange in a Multielectron Quantum Dot, Phys. Rev. Lett. 119, 227701 (2017).
  * Malinowski  _et al._ [2018] F. K. Malinowski, F. Martins, T. B. Smith, S. D. Bartlett, A. C. Doherty, P. D. Nissen, S. Fallahi, G. C. Gardner, M. J. Manfra, C. M. Marcus, and F. Kuemmeth, Spin of a multielectron quantum dot and its interaction with a neighboring electron, Phys. Rev. X 8, 011045 (2018).
  * Wang  _et al._ [2023] Z. Wang, M. Feng, S. Serrano, W. Gilbert, R. C. C. Leon, T. Tanttu, P. Mai, D. Liang, J. Y. Huang, Y. Su, W. H. Lim, F. E. Hudson, C. C. Escott, A. Morello, C. H. Yang, A. S. Dzurak, A. Saraiva, and A. Laucht, Jellybean quantum dots in silicon for qubit coupling and on-chip quantum chemistry, Advanced Materials 35, 2208557 (2023).
  * Patomäki  _et al._ [2024] S. M. Patomäki, J. Williams, F. Berritta, C. Lainé, M. A. Fogarty, R. C. C. Leon, J. Jussot, S. Kubicek, A. Chatterjee, B. Govoreanu, F. Kuemmeth, J. J. L. Morton, and M. F. Gonzalez-Zalba, Elongated quantum dot as a distributed charge sensor, Phys. Rev. Appl. 21, 054042 (2024).
  * Note [1] In the case of a confinement with circular symmetry, the WMs are referred to as rotating WMs as well.
  * Yannouleas and Landman [1999] C. Yannouleas and U. Landman, Spontaneous Symmetry Breaking in Single and Molecular Quantum Dots, Phys. Rev. Lett. 82, 5325–5328 (1999).
  * Yannouleas and Landman [2000] C. Yannouleas and U. Landman, Collective and Independent-Particle Motion in Two-Electron Artificial Atoms, Phys. Rev. Lett. 85, 1726–1729 (2000).
  * Yannouleas and Landman [2002] C. Yannouleas and U. Landman, Strongly correlated wavefunctions for artificial atoms and molecules, Journal of Physics: Condensed Matter 14, L591–L598 (2002).
  * Harju  _et al._ [2002] A. Harju, S. Siljamäki, and R. M. Nieminen, Wigner molecules in quantum dots: A quantum Monte Carlo study, Phys. Rev. B 65, 075309 (2002).
  * Mikhailov [2002] S. A. Mikhailov, Two ground-state modifications of quantum-dot beryllium, Phys. Rev. B 66, 153313 (2002).
  * Yannouleas and Landman [2003] C. Yannouleas and U. Landman, Two-dimensional quantum dots in high magnetic fields: Rotating-electron-molecule versus composite-fermion approach, Phys. Rev. B 68, 035326 (2003).
  * Tavernier  _et al._ [2003] M. B. Tavernier, E. Anisimovas, F. M. Peeters, B. Szafran, J. Adamowski, and S. Bednarek, Four-electron quantum dot in a magnetic field, Phys. Rev. B 68, 205305 (2003).
  * Yannouleas and Landman [2007] C. Yannouleas and U. Landman, Symmetry breaking and quantum correlations in finite systems: Studies of quantum dots and ultracold Bose gases and related nuclear and chemical methods, Reports on Progress in Physics 70, 2067–2148 (2007).
  * Yannouleas and Landman [2011] C. Yannouleas and U. Landman, Unified microscopic approach to the interplay of pinned-Wigner-solid and liquid behavior of the lowest Landau-level states in the neighborhood of $\nu=1/3$, Phys. Rev. B 84, 165327 (2011).
  * Ellenberger  _et al._ [2006] C. Ellenberger, T. Ihn, C. Yannouleas, U. Landman, K. Ensslin, D. Driscoll, and A. C. Gossard, Excitation Spectrum of Two Correlated Electrons in a Lateral Quantum Dot with Negligible Zeeman Splitting, Phys. Rev. Lett. 96, 126806 (2006).
  * Ihn  _et al._ [2007] T. Ihn, C. Ellenberger, K. Ensslin, C. Yannouleas, U. Landman, D. C. Driscoll, and A. C. Gossard, Quantum dots based on parabolic quantum wells: Importance of electronic correlations, International Journal of Modern Physics B 21, 1316–1325 (2007).
  * Jang  _et al._ [2021] W. Jang, M.-K. Cho, H. Jang, J. Kim, J. Park, G. Kim, B. Kang, H. Jung, V. Umansky, and D. Kim, Single-Shot Readout of a Driven Hybrid Qubit in a GaAs Double Quantum Dot, Nano Letters 21, 4999–5005 (2021).
  * Yannouleas and Landman [2022a] C. Yannouleas and U. Landman, Wigner molecules and hybrid qubits, J. Phys.: Condens. Matter (Letter) 34, 21LT01 (2022a).
  * Yannouleas and Landman [2022b] C. Yannouleas and U. Landman, Molecular formations and spectra due to electron correlations in three-electron hybrid double-well qubits, Phys. Rev. B 105, 205302 (2022b).
  * Corrigan  _et al._ [2021] J. Corrigan, J. P. Dodson, H. E. Ercan, J. C. Abadillo-Uriel, B. Thorgrimsson, T. J. Knapp, N. Holman, T. McJunkin, S. F. Neyens, E. R. MacQuarrie, R. H. Foote, L. F. Edge, M. Friesen, S. N. Coppersmith, and M. A. Eriksson, Coherent control and spectroscopy of a semiconductor quantum dot Wigner molecule, Phys. Rev. Lett. 127, 127701 (2021).
  * Yannouleas and Landman [2022c] C. Yannouleas and U. Landman, Valleytronic full configuration-interaction approach: Application to the excitation spectra of Si double-dot qubits, Phys. Rev. B 106, 195306 (2022c).
  * Pecker  _et al._ [2013] S. Pecker, F. Kuemmeth, A. Secchi, M. Rontani, D. C. Ralph, P. L. McEuen, and S. Ilani, Observation and spectroscopy of a two-electron Wigner molecule in an ultraclean carbon nanotube, Nature Physics 9, 576–581 (2013).
  * Yannouleas and Landman [2023] C. Yannouleas and U. Landman, Quantum Wigner molecules in moiré materials, Phys. Rev. B 108, L121411 (2023).
  * Yannouleas and Landman [2024] C. Yannouleas and U. Landman, Wigner-molecule supercrystal in transition metal dichalcogenide moiré superlattices: Lessons from the bottom-up approach, Phys. Rev. B 109, L121302 (2024).
  * [28] H. Li, Z. Xiang, A. P. Reddy, T. Devakul, R. Sailus, R. Banerjee, T. Taniguchi, K. Watanabe, S. Tongay, A. Zettl, L. Fu, M. F. Crommie, and F. Wang, Wigner molecular crystals from multi-electron moiré artificial atoms, arXiv:2312.07607 .
  * Löwdin [1955] P.-O. Löwdin, Quantum theory of many-particle systems. I. Physical interpretations by means of density matrices, natural spin-orbitals, and convergence problems in the method of configurational interaction, Phys. Rev. 97, 1474–1489 (1955).
  * Shavitt [1998] I. Shavitt, The history and evolution of configuration interaction, Molecular Physics 94, 3–17 (1998).
  * Szabo and Ostlund [1989] A. Szabo and N. S. Ostlund,  _Modern Quantum Chemistry_ (McGraw-Hill, New York, 1989).
  * Piacente  _et al._ [2004] G. Piacente, I. V. Schweigert, J. J. Betouras, and F. M. Peeters, Generic properties of a quasi-one-dimensional classical Wigner crystal, Phys. Rev. B 69, 045324 (2004).
  * Meyer and Matveev [2008] J. S. Meyer and K. A. Matveev, Wigner crystal physics in quantum wires, Journal of Physics: Condensed Matter 21, 023203 (2008).
  * Ho  _et al._ [2018] S.-C. Ho, H.-J. Chang, C.-H. Chang, S.-T. Lo, G. Creeth, S. Kumar, I. Farrer, D. Ritchie, J. Griffiths, G. Jones, M. Pepper, and T.-M. Chen, Imaging the zigzag Wigner crystal in confinement-tunable quantum wires, Phys. Rev. Lett. 121, 106801 (2018).
  * Kumar and Pepper [2021] S. Kumar and M. Pepper, Chapter Two - Advances in interaction effects in the quasi one-dimensional electron gas, in  _Semiconductor Nanodevices_, Frontiers of Nanoscience, Vol. 20, edited by D. A. Ritchie (Elsevier, 2021) pp. 7–29.
  * Sheikh  _et al._ [2021] J. A. Sheikh, J. Dobaczewski, P. Ring, L. M. Robledo, and C. Yannouleas, Symmetry restoration in mean-field approaches, Journal of Physics G: Nuclear and Particle Physics 48, 123001 (2021).
  * Zwanenburg  _et al._ [2013b] F. A. Zwanenburg, A. S. Dzurak, A. Morello, M. Y. Simmons, L. C. L. Hollenberg, G. Klimeck, S. Rogge, S. N. Coppersmith, and M. A. Eriksson, Silicon quantum electronics, Rev. Mod. Phys. 85, 961–1019 (2013b).
  * Note [2] These confinement-induced single-particle energy states are also referred to as “orbitals”, see, e.g., the expressions “atomic orbitals”, “space orbitals”, and “spin-orbitals” in chemistry [31].
  * Note [3] This is an apparent generalization of the term spin-orbital used in chemistry and molecular physics.
  * Note [4] This behavior contrasts with that in harmonic confinements along the lateral $x$-direction, where the zig-zag chain starts forming at the center of the linear chain, a fact well known from classical calculations in the literature of trapped heavy ions; see, e.g., Refs. [46, 44, 45].
  * Pacati and Boffi [1970] F. D. Pacati and S. Boffi, Slater determinants, parity projection, and Hartree-Fock calculations, Phys. Rev. C 2, 1205–1210 (1970).
  * Zhang  _et al._ [2023] J. Zhang, B. T. Chow, S. Ejtemaee, and P. C. Haljan, Spectroscopic characterization of the quantum linear-zigzag transition in trapped ions, npj Quantum Information 9, 68 (2023).
  * Li  _et al._ [2024] H. Li, Z. Xiang, T. Wang, M. H. Naik, W. Kim, J. Nie, S. Li, Z. Ge, Z. He, Y. Ou, R. Banerjee, T. Taniguchi, K. Watanabe, S. Tongay, A. Zettl, S. G. Louie, M. P. Zaletel, M. F. Crommie, and F. Wang, Imaging tunable Luttinger liquid systems in van der Waals heterostructures (2024), arXiv:2404.16344 [cond-mat.str-el] .
  * Ejtemaee [2015] S. Ejtemaee,  _Dynamics of Trapped Ions Near the Linear–Zigzag Structural Phase Transition_ , Phd thesis, Simon Fraser University (2015), available at https://summit.sfu.ca/item/16248.
  * Yan  _et al._ [2016] L. L. Yan, W. Wan, L. Chen, F. Zhou, S. J. Gong, X. Tong, and M. Feng, Exploring structural phase transitions of ion crystals, Scientific Reports 6, 21547 (2016).
  * Schiffer [1993] J. P. Schiffer, Phase transitions in anisotropically confined ionic crystals, Phys. Rev. Lett. 70, 818–821 (1993).