# Classical modelling of a lossy Gaussian bosonic sampler

M.V. Umanskii  Department of Physics, Lomonosov Moscow State University,
Leninskie gory 1, Moscow 119991, Russia  A.N. Rubtsov  Department of Physics,
Lomonosov Moscow State University, Leninskie gory 1, Moscow 119991, Russia
Russian Quantum Center, Bolshoy Bulvar 30, bld. 1, Skolkovo, Moscow 121205,
Russia

###### Abstract

Gaussian boson sampling (GBS) is considered a candidate problem for
demonstrating quantum advantage. We propose an algorithm for approximate
classical simulation of a lossy GBS instance. The algorithm relies on the
Taylor series expansion, and increasing the number of terms of the expansion
that are used in the calculation yields greater accuracy. The complexity of
the algorithm is polynomial in the number of modes given the number of terms
is fixed. We describe conditions for the input state squeezing parameter and
loss level that provide the best efficiency for this algorithm (by efficient
we mean that the Taylor series converges quickly). In recent experiments that
claim to have demonstrated quantum advantage, these conditions are satisfied;
thus, this algorithm can be used to classically simulate these experiments.

##  1 Introduction

Quantum computers are computational devices which operate using phenomena
described by quantum mechanics. Therefore, they can carry out the operations
which are not available for classical computers. The ability of a quantum
computer to solve a specific task faster that any classical computer is
usually referred to as quantum advantage. Although quantum algorithms that
provide exponential speedup over classical ones are known, they are hard to
implement in practice. Examples of such algorithms include Shor’s algorithm of
factoring integers[1], that works in polynomial time, whereas all classical
algorithms require exponential time. Modern quantum computers are far from
experimentally demonstrating quantum advantage on basic problems like integer
factorization.

Boson sampling[2] is a problem that was proposed as a good candidate for
demonstrating quantum advantage due to its nature. A boson sampler is a
linear-optical device that consists of non-classical sources of
indistinguishable photons, a multichannel interferometer mixing photons of
different sources, and photon detectors at the output channels of the
interferometer. In the original proposal, the indistinguishable photons were
prepared in Fock states. The problem then is to calculate the photon
statistics after the interferometer, given an input state and the
interferometer matrix. The relevant parameters are the number of modes $N$ and
the total number of photons injected in the interferometer $M$. Experimentally
it corresponds to performing multiple measurements of the photon counts at the
outputs of such a device[3].

Due to the technological complexity of generating Fock states, several
variants of the original boson sampling problem have been proposed. They aim
at improving the photon generation efficiency and increasing the scale of
implementations. One such example is the Scattershot boson sampling, which
uses many parametric down-conversion sources to improve the single photon
generation rate. It has been implemented experimentally using a 13-mode
integrated photonic chip and six PDC photon sources[4].

Another variant is the Gaussian boson sampling[5][6], in which Gaussian states
are injected into the interferometer instead of Fock states. Gaussian input
states can be generated using PDC sources, and it allows the non-classical
input states to be prepared deterministically. In this variant, the relative
input photon phases can affect the sampling distribution. Experiments were
carried out with $N=12$[7], $N=100$[8] and $N=144$[9]. The latter
implementations used PPKTP crystals as PDC sources and employs an active phase
locking mechanism to ensure a coherent superposition.

Any experimental set-up, of course, differs from the idealized model
considered in theoretical modelling. Bosonic samplers suffer from two
fundamental types of imperfections. First, the parameters of a real device,
such as the reflection coefficients of the beam splitters and the phase
rotations, are never known exactly. A small change in the interferometer
parameters can affect the sampling statistics drastically, so that modelling
of an ideal device no longer makes much sense. Another type of imperfections
is photon losses. These losses happen because of imperfections in photon
preparation, absorption inside the interferometer and imperfect detectors and
coupling.

There are different ways of modelling losses, for example by introducing extra
beam splitters[10] or replacing the interferometer matrix by a combination of
lossless linear optics transformations and the diagonal matrix that contains
transmission coefficients[11]. In the algorithm described in this paper we
will assume that losses occur on the inputs of the interferometer, and we will
describe the exact way that we model them.

Imperfections in middle-sized systems make them, in general, easier to emulate
with classical computers[12]. It was shown[13] that with the increase of
losses in a system the complexity of the task decreases. When the number of
photons $M^{\prime}$ that arrive at the outputs is less than $\sqrt{M}$, the
problem of boson sampling can be efficiently solved using classical computers.
On the other hand, if the losses are low, the problem remains hard for
classical computers[14].

In this paper we propose a classical algorithm for calculating probabilities
of output states in a GBS problem. The algorithm uses Taylor series expansion,
and it converges faster depending on the parameters of the problem, namely the
amount of losses in the system and the squeezing parameter of the input
states. The higher the losses in the system, the less orders of the series are
needed to approximate the probability of observing a given output state.

The work by Oh et al.[15] used the following approach to simulating GBS: the
covariance matrix of the output Gaussian state was decomposed into ”quantum”
and ”classical” parts, and the ”quantum” part was simulated using matrix
product states and the ”classical” part was simulated by random displacement.
Thus, when the photon loss rate is high, computational complexity of this
algorithm is reduced.

The algorithm that we propose in this paper uses some similar ideas, namely
the zeroth order of the Taylor series may be considered the ”classical” part
that is computed quite easily, while the remaining terms are the ”quantum”
part that is more computationally complex. The contribution of this ”quantum”
part is smaller when the losses in the system are high; thus, our algorithm
also has optimal conditions that depend on the magnitude of losses. We also
analyze some recent GBS implementations to compare the conditions in those
experiments with the optimal conditions for our algorithm.

##  2 Problem specification

Let us first consider a lossless linear-optics interferometer with a
transmission matrix $U$:

$\hat{a}_{i}^{\dagger}=\sum_{j}U_{ij}\hat{d}_{j}^{\dagger},~{}~{}~{}\hat{a}_{i}=\sum_{j}U_{ij}^{*}\hat{d}_{j}$
(1)

where creation operators acting on the $i$-th input and output modes are
denoted $a_{i}^{\dagger}$ and $d_{i}^{\dagger}$. Suppose the input modes are
injected with single-mode squeezed states:

$\ket{\psi}=e^{\sum_{i}\frac{\alpha_{i}}{2}(\hat{a}_{i}^{\dagger})^{2}}\ket{0}.$
(2)

The goal is to calculate the probability of detecting $n_{1}$ photons in the
first output mode, $n_{2}$ photons in the second output mode and so on. This
probability can be calculated in the following way:

$Tr\left\\{\hat{\rho}_{out}\hat{\vec{n}}\right\\}=Tr\left\\{\hat{\rho}_{out}\bigotimes_{i}\ket{n_{i}}\bra{n_{i}}\right\\},$
(3)

where $\hat{\rho}_{out}$ is the density matrix of the output state.

###  2.1 Modelling losses

In real-life bosonic samplers there will always be losses. Here we will model
them by substituting

$a_{i}^{\dagger}\longrightarrow ca_{i}^{\dagger}+sb_{i}^{\dagger},$ (4)

where $b_{i}^{\dagger}$ acts on a mode that we cannot observe, and
$c^{2}+s^{2}=1$, $c,s\in\mathbb{R}$. Now, the goal is to compute the same
probability ($Tr\left\\{\hat{\rho_{out}}\hat{\vec{n}}\right\\}$), but taking
losses into account. Input state will now be

$\ket{\psi^{\prime}}=e^{\sum_{i}\frac{\alpha_{i}}{2}(c\hat{a}_{i}^{\dagger}+s\hat{d}_{i}^{\dagger})^{2}}\ket{0_{a}0_{b}},$
(5)

and we now take partial trace over all loss modes when calculating the density
matrix:

$\hat{\rho}=Tr_{b}\left\\{e^{\sum_{i}\frac{\alpha_{i}}{2}(c\hat{a}_{i}^{\dagger}+s\hat{d}_{i}^{\dagger})^{2}}\ket{0_{a}0_{b}}\bra{0_{a}0_{b}}e^{\sum_{i}\frac{\alpha_{i}}{2}(c\hat{a}_{i}^{\dagger}+s\hat{d}_{i}^{\dagger})^{2}}\right\\}.$
(6)

##  3 Algorithm derivation

Let us consider a single mode:

$\ket{\psi^{\prime}}=e^{\frac{\alpha}{2}(c\hat{a}^{\dagger}+s\hat{b}^{\dagger})^{2}}\ket{0_{a}0_{b}},$
(7)
$\hat{\rho}=Tr_{b}\left\\{e^{\frac{\alpha}{2}(c\hat{a}^{\dagger}+s\hat{b}^{\dagger})^{2}}\ket{0_{a}0_{b}}\bra{0_{a}0_{b}}e^{\frac{\alpha}{2}(c\hat{a}+s\hat{b})^{2}}\right\\}.$
(8)

###  3.1 Calculating partial trace

We start by applying the Hubbard–Stratonovich transformation[16][17]

$e^{\frac{\hat{A}^{2}}{2}}=\frac{1}{\sqrt{2\pi}}\int
e^{\xi\hat{A}-\frac{\xi^{2}}{2}}d\xi$ (9)

to both exponents in the density matrix operator. This gives us the following:

$\hat{\rho}=\frac{1}{2\pi}\int
Tr_{b}\left\\{e^{\xi\sqrt{\alpha}(c\hat{a}^{\dagger}+s\hat{b}^{\dagger})}\ket{0_{a}0_{b}}\bra{0_{a}0_{b}}e^{\tilde{\xi}\sqrt{\alpha}(c\hat{a}+s\hat{b})}\right\\}e^{-\frac{\xi^{2}+\tilde{\xi}^{2}}{2}}d\xi
d\tilde{\xi}=$ $=\frac{1}{2\pi\alpha}\int
Tr_{b}\left\\{e^{\xi\sqrt{\alpha}(c\hat{a}^{\dagger}+s\hat{b}^{\dagger})}\ket{0_{a}0_{b}}\bra{0_{a}0_{b}}e^{\tilde{\xi}\sqrt{\alpha}(c\hat{a}+s\hat{b})}\right\\}e^{-\frac{(\xi\sqrt{\alpha})^{2}+(\tilde{\xi}\sqrt{\alpha})^{2}}{2\alpha}}d(\xi\sqrt{\alpha})d(\tilde{\xi}\sqrt{\alpha}).$

Let us redefine $\xi\sqrt{\alpha}\longrightarrow\xi$,
$\tilde{\xi}\sqrt{\alpha}\longrightarrow\tilde{\xi}$ for convenience:

$\hat{\rho}=\frac{1}{2\pi\alpha}\int
Tr_{b}\left\\{e^{\xi(c\hat{a}^{\dagger}+s\hat{b}^{\dagger})}\ket{0_{a}0_{b}}\bra{0_{a}0_{b}}e^{\tilde{\xi}(c\hat{a}+s\hat{b})}\right\\}e^{-\frac{\xi^{2}+\tilde{\xi}^{2}}{2\alpha}}d\xi
d\tilde{\xi}.$ (10)

We can now calculate the partial trace over loss modes:

$Tr_{b}\left\\{e^{\xi(c\hat{a}^{\dagger}+s\hat{b}^{\dagger})}\ket{0_{a}0_{b}}\bra{0_{a}0_{b}}e^{\tilde{\xi}(c\hat{a}+s\hat{b})}\right\\}=$
$=e^{\xi c\hat{a}^{\dagger}}\ket{0_{a}}\bra{0_{a}}e^{\tilde{\xi}c\hat{a}}\cdot
Tr\left\\{e^{\xi
s\hat{b}^{\dagger}}\ket{0_{b}}\bra{0_{b}}e^{\tilde{\xi}s\hat{b}}\right\\}=$
$=e^{\xi
c\hat{a}^{\dagger}}\ket{0_{a}}\bra{0_{a}}e^{\tilde{\xi}c\hat{a}}\cdot\bra{0_{b}}e^{\tilde{\xi}s\hat{b}}e^{\xi
s\hat{b}^{\dagger}}\ket{0_{b}}.$

The following expression can be simplified:

$\bra{0_{b}}e^{\tilde{\xi}s\hat{b}}e^{\xi s\hat{b}^{\dagger}}\ket{0_{b}}=$
$=\bra{0_{b}}(1+\tilde{\xi}s\hat{b}+\frac{1}{2}(\xi s\hat{b})^{2}+...)(1+\xi
s\hat{b^{\dagger}}+\frac{1}{2}(\xi s\hat{b^{\dagger}})^{2}+...)\ket{0_{b}}=$
$=\left(\bra{0_{b}}+\tilde{\xi}s\bra{1_{b}}+\frac{1}{\sqrt{2}}(\tilde{\xi}s)^{2}\bra{2_{b}}+...\right)\left(\ket{0_{b}}+\xi
s\ket{1_{b}}+\frac{1}{\sqrt{2}}(\xi s)^{2}\ket{2_{b}}+...\right)=$
$=1+\xi\tilde{\xi}s^{2}+\frac{1}{2}(\xi\tilde{\xi}s^{2})^{2}+...=e^{\xi\tilde{\xi}s^{2}}.$

The density matrix now can be written in the following way:

$\hat{\rho}=\frac{1}{2\pi\alpha}\int e^{\xi
c\hat{a}^{\dagger}}\ket{0}\bra{0}e^{\tilde{\xi}c\hat{a}}\cdot
e^{-\frac{\xi^{2}+\tilde{\xi}^{2}}{2\alpha}+\xi\tilde{\xi}s^{2}}d\xi
d\tilde{\xi}.$ (11)

###  3.2 Switching between probability density functions

We can view this integral as taking an expected value over a $2$-dimensional
normal distribution. $\xi$ and $\tilde{\xi}$ then become normally distributed
random variables with mean vector equal to zero. Their covariance matrix has
the following form:

$\Sigma=\begin{pmatrix}1/\alpha&-s^{2}\\\ -s^{2}&1/\alpha\\\
\end{pmatrix}^{-1}=\frac{1}{1/\alpha^{2}-s^{4}}\begin{pmatrix}1/\alpha&s^{2}\\\
s^{2}&1/\alpha\\\ \end{pmatrix}.$ (12)

Then we can write

$\hat{\rho}=\frac{(det\Sigma)^{1/2}}{\alpha}\frac{1}{2\pi(det\Sigma)^{1/2}}\int
e^{\xi
c\hat{a}^{\dagger}}\ket{0}\bra{0}e^{\tilde{\xi}c\hat{a}}e^{-\frac{\xi^{2}+\tilde{\xi}^{2}}{2\alpha}+\xi\tilde{\xi}s^{2}}d\xi
d\tilde{\xi}=$ (13)
$=\frac{(det\Sigma)^{1/2}}{\alpha}\cdot\mathbb{E}_{\mathbb{N}(0,\Sigma)}\left[e^{\xi
c\hat{a}^{\dagger}}\ket{0}\bra{0}e^{\tilde{\xi}c\hat{a}}\right],$ (14)

where $\mathbb{E}_{\mathbb{N}(0,\Sigma)}$ denotes averaging over the
$2$-dimensional normal distribution $\mathbb{N}(0,\Sigma)$.

The expression $e^{\xi
c\hat{a}^{\dagger}}\ket{0}\bra{0}e^{\tilde{\xi}c\hat{a}}$ is troublesome to
calculate, since there are two different variables $\xi$ and $\tilde{\xi}$. We
want to arrive somehow at an expression with only one such variable, i.e.
$e^{\xi c\hat{a}^{\dagger}}\ket{0}\bra{0}e^{\xi c\hat{a}}$, which we will
denote $\hat{\nu}(\xi c)$.

We now will choose normally distributed random variables
$\xi_{0},\chi,\tilde{\chi}\in\mathbb{R}$ such that
$\xi=\xi_{0}+\chi,~{}\tilde{\xi}=\xi_{0}+\tilde{\chi}$ and the distributions
over $\xi,\tilde{\xi}$ and $\xi_{0},\chi,\tilde{\chi}$ have the same moments:

$\begin{cases}\overline{\xi^{2}}=\overline{(\xi_{0}+\chi)^{2}}=\overline{\xi_{0}^{2}}+2\overline{\xi_{0}\chi}+\overline{\chi^{2}},\\\
\overline{\tilde{\xi}^{2}}=\overline{(\xi_{0}+\tilde{\chi})^{2}}=\overline{\xi_{0}^{2}}+2\overline{\xi_{0}\tilde{\chi}}+\overline{\tilde{\chi}^{2}},\\\
\overline{\xi\tilde{\xi}}=\overline{(\xi_{0}+\chi)(\xi_{0}+\tilde{\chi})}=\overline{\xi_{0}^{2}}+\overline{\xi_{0}\chi}+\overline{\xi_{0}\tilde{\chi}}+\overline{\chi\tilde{\chi}}.\end{cases}$
(15)

We have some freedom in choosing these variables; we will set
$\overline{\xi_{0}\chi}=\overline{\xi_{0}\tilde{\chi}}=0$ so that
$\xi_{0}\perp\\!\\!\\!\perp\chi$ and $\xi_{0}\perp\\!\\!\\!\perp\tilde{\chi}$.
Then the covariance matrix $\Gamma$ of $\xi_{0},\chi,\tilde{\chi}$ will be
determined by one parameter $h=\overline{\chi\tilde{\chi}}$:

$\begin{cases}\overline{\xi^{2}}=\overline{\xi_{0}^{2}}+\overline{\chi^{2}},\\\
\overline{\tilde{\xi}^{2}}=\overline{\xi_{0}^{2}}+\overline{\tilde{\chi}^{2}},\\\
\overline{\xi\tilde{\xi}}=\overline{\xi_{0}^{2}}+h.\end{cases}$ (16)

$\begin{cases}\overline{\xi_{0}^{2}}=\overline{\xi\tilde{\xi}}-h=\frac{s^{2}}{1/\alpha^{2}-s^{4}}-h,\\\
\overline{\chi^{2}}=\overline{\tilde{\chi}^{2}}=\overline{\xi^{2}}-\overline{\xi\tilde{\xi}}+h=\frac{1/\alpha-s^{2}}{1/\alpha^{2}-s^{4}}+h=\frac{1}{1/\alpha+s^{2}}+h.\end{cases}$
(17)

Note that $-\frac{1}{1/\alpha+s^{2}}\leq
h\leq\frac{s^{2}}{1/\alpha^{2}-s^{4}}$. We will later find an optimal way to
choose $h$. The density matrix in terms of the new variables
$\xi_{0},\chi,\tilde{\chi}\in\mathbb{N}(0,\Gamma)$ is

$\hat{\rho}=\frac{(det\Sigma)^{1/2}}{\alpha}\cdot\mathbb{E}_{\mathbb{N}(0,\Gamma)}\left[e^{(\xi_{0}+\chi)c\hat{a}^{\dagger}}\ket{0}\bra{0}e^{(\xi_{0}+\tilde{\chi})c\hat{a}}\right].$
(18)

###  3.3 Taylor series expansion

We now consider the Taylor series of the expression
$e^{(\xi_{0}+\chi)c\hat{a}^{\dagger}}\ket{0}\bra{0}e^{(\xi_{0}+\tilde{\chi})c\hat{a}}$,
leaving only $\xi_{0}$ in the exponent:

$e^{(\xi_{0}+\chi)c\hat{a}^{\dagger}}\ket{0}\bra{0}e^{(\xi_{0}+\tilde{\chi})c\hat{a}}=e^{\chi
c\hat{a}^{\dagger}}e^{\xi_{0}c\hat{a}^{\dagger}}\ket{0}\bra{0}e^{\xi_{0}c\hat{a}}e^{\tilde{\chi}c\hat{a}}=$
$=e^{\chi
c\hat{a}^{\dagger}}\hat{\nu}(\xi_{0}c)e^{\tilde{\chi}c\hat{a}}=\left(1+\chi
c\hat{a}^{\dagger}+\frac{(\chi
c\hat{a}^{\dagger})^{2}}{2}+...\right)\hat{\nu}(\xi_{0}c)\left(1+\tilde{\chi}c\hat{a}+\frac{(\tilde{\chi}c\hat{a})^{2}}{2}+...\right).$

Each term in the expression will be proportional to

$\chi^{n}\tilde{\chi}^{m}\cdot(\hat{a}^{\dagger})^{n}\hat{\nu}(\xi_{0}c)\hat{a}^{m},$

and since $\xi_{0}\perp\\!\\!\\!\perp\chi$ and
$\xi_{0}\perp\\!\\!\\!\perp\tilde{\chi}$, the integral over
$\xi_{0},\chi,\tilde{\chi}$ can be written as a product of integrals over
$\xi_{0}$ and $\chi,\tilde{\chi}$. The latter can be taken analytically:

$\int\chi^{n}\tilde{\chi}^{m}d\chi
d\tilde{\chi}\propto\mathbb{E}_{\mathbb{N}(0,\Gamma)}\left[\chi^{n}\tilde{\chi}^{m}\right],$
(19)

which can be calculated using Wick’s probability theorem.

###  3.4 Choosing $\Gamma$

The idea consists in minimizing the ”perturbation parameter” so that each
subsequent order of the Taylor series expansion has less impact on the
expression. Since higher orders of the expansion contain higher powers of
$c^{2}$ and higher moments
$\mathbb{E}_{\mathbb{N}(0,\Gamma)}\left[\chi^{n}\tilde{\chi}^{m}\right]$, and
these moments can be calculated via second moments
$\overline{\chi^{2}}=\overline{\tilde{\chi}^{2}}$ and
$\overline{\chi\tilde{\chi}}=h$, the role of the ”perturbation parameter” is
played by
$\varepsilon=c^{2}\cdot\max(\overline{\chi^{2}},|\overline{\chi\tilde{\chi}}|)$.

Let us consider the conditions that must be satisfied by $h$. Firstly, $h$
must satisfy $-\frac{1}{1/\alpha+s^{2}}\leq
h\leq\frac{s^{2}}{1/\alpha^{2}-s^{4}}$, because $\overline{\xi_{0}^{2}}\geq 0$
and $\overline{\chi^{2}}\geq 0$. Secondly, since $\Gamma$ is a covariance
matrix, its eigenvalues must be non-negative. The eigenvalues of $\Gamma$ are
$\overline{\xi_{0}^{2}}$, $\overline{\chi^{2}}-h$ and $\overline{\chi^{2}}+h$.
Thus, $h$ needs to satisfy

$\overline{\chi^{2}}+h\geq 0\iff\frac{1}{1/\alpha+s^{2}}+2h\geq 0\iff
h\geq-\frac{1}{2}\frac{1}{1/\alpha+s^{2}}.$ (20)

The minimum of $\max(\overline{\chi^{2}},|h|)$ is realised when
$h=-\overline{\chi^{2}}=-\frac{1}{2}\frac{1}{1/\alpha+s^{2}}$.

###  3.5 Multimode case

Let’s apply the steps described above to the case of $N$ modes. We start with
an input state

$\ket{\psi^{\prime(N)}}=\prod_{i=1}^{N}e^{\frac{\alpha}{2}(c\hat{a}_{i}^{\dagger}+s\hat{b}_{i}^{\dagger})^{2}}\ket{0_{a}0_{b}}.$

We construct a density matrix and take the partial trace over loss modes:

$\hat{\rho}=Tr_{b}\left\\{e^{\sum_{i}\frac{\alpha}{2}(c\hat{a}_{i}^{\dagger}+s\hat{b}_{i}^{\dagger})^{2}}\ket{0_{a}0_{b}}\bra{0_{a}0_{b}}e^{\sum_{i}\frac{\alpha}{2}(c\hat{a}_{i}+s\hat{b}_{i})^{2}}\right\\}.$

We apply the Hubbard–Stratonovich transformation $2N$ times, resulting in an
integral over $\prod_{i=1}^{N}d\xi_{i}d\tilde{\xi}_{i}$:

$\hat{\rho}=\frac{1}{(2\pi)^{N}}\int
Tr_{b}\left\\{e^{\sum_{i}\xi_{i}\sqrt{\alpha}(c\hat{a}_{i}^{\dagger}+s\hat{b}_{i}^{\dagger})}\ket{0_{a}0_{b}}\bra{0_{a}0_{b}}e^{\sum_{i}\tilde{\xi}_{i}\sqrt{\alpha}(c\hat{a}_{i}+s\hat{b}_{i})}\right\\}e^{-\sum_{i}\frac{\xi_{i}^{2}+\tilde{\xi}_{i}^{2}}{2}}\prod_{i}d\xi_{i}d\tilde{\xi}_{i}=$
$=\frac{1}{(2\pi\alpha)^{N}}\int
Tr_{b}\left\\{e^{\sum_{i}\xi_{i}\sqrt{\alpha}(c\hat{a}_{i}^{\dagger}+s\hat{b}_{i}^{\dagger})}\ket{0_{a}0_{b}}\bra{0_{a}0_{b}}e^{\sum_{i}\tilde{\xi}_{i}\sqrt{\alpha}(c\hat{a}_{i}+s\hat{b}_{i})}\right\\}\cdot$
$\cdot
e^{-\sum_{i}\frac{(\xi_{i}\sqrt{\alpha})^{2}+(\tilde{\xi}_{i}\sqrt{\alpha})^{2}}{2\alpha}}\prod_{i}d(\xi_{i}\sqrt{\alpha})d(\tilde{\xi}_{i}\sqrt{\alpha}).$

Again, we redefine $\xi_{i}\sqrt{\alpha}\longrightarrow\xi_{i}$,
$\tilde{\xi}_{i}\sqrt{\alpha}\longrightarrow\tilde{\xi}_{i}$:

$\hat{\rho}=\frac{1}{(2\pi\alpha)^{N}}\int
Tr_{b}\left\\{e^{\sum_{i}\xi_{i}(c\hat{a}_{i}^{\dagger}+s\hat{b}_{i}^{\dagger})}\ket{0_{a}0_{b}}\bra{0_{a}0_{b}}e^{\sum_{i}\tilde{\xi}_{i}(c\hat{a}_{i}+s\hat{b}_{i})}\right\\}e^{-\sum_{i}\frac{\xi_{i}^{2}+\tilde{\xi}_{i}^{2}}{2\alpha_{i}}}\prod_{i}d\xi_{i}d\tilde{\xi}_{i}.$

We compute partial trace over loss modes:

$\hat{\rho}=\frac{1}{(2\pi\alpha)^{N}}\int
e^{\sum_{i}\xi_{i}c\hat{a}_{i}^{\dagger}}\ket{0}\bra{0}e^{\sum_{i}\tilde{\xi}_{i}c\hat{a}_{i}}e^{-\sum_{i}\frac{\xi_{i}^{2}+\tilde{\xi}_{i}^{2}}{2\alpha}+\xi_{i}\tilde{\xi}_{i}s^{2}}\prod_{i}d\xi_{i}d\tilde{\xi}_{i}.$

This expression now can be considered as taking an expected value over a
$2N$-dimensional normal distribution, where all variable pairs
$\xi_{i},\tilde{\xi}_{i}$ are independent. Every variable pair
$\xi_{i},\tilde{\xi}_{i}$ has covariance matrix $\Sigma$, and we can write
this expression in the following way:

$\hat{\rho}=\frac{(det\Sigma)^{N/2}}{\alpha^{N}}\cdot\mathbb{E}_{\prod_{i}\mathbb{N}(0,\Sigma)}\left[e^{\sum_{i}\xi_{i}c\hat{a}_{i}^{\dagger}}\ket{0}\bra{0}e^{\sum_{i}\tilde{\xi_{i}}c\hat{a}_{i}}\right].$
(21)

For each variable pair $\xi_{i},\tilde{\xi}_{i}$ we now choose
$\xi_{0i},\chi_{i},\tilde{\chi}_{i}$ in a way that is described above. Then,

$\hat{\rho}=\frac{(det\Sigma)^{N/2}}{\alpha^{N}}\cdot\mathbb{E}_{\prod_{i}\mathbb{N}(0,\Gamma)}\left[e^{\sum_{i}(\xi_{0i}+\chi_{i})c\hat{a}_{i}^{\dagger}}\ket{0}\bra{0}e^{\sum_{i}(\xi_{0i}+\tilde{\chi}_{i})c\hat{a}_{i}}\right].$
(22)

We now consider the Taylor series expansion (up to the second order) of the
expression in the square brackets, which we will denote $\hat{\mu}$:

$\hat{\mu}=e^{\sum_{i}\chi_{i}c\hat{a}_{i}^{\dagger}}e^{\sum_{i}\xi_{0i}c\hat{a}_{i}^{\dagger}}\ket{0}\bra{0}e^{\sum_{i}\xi_{0i}c\hat{a}_{i}}e^{\sum_{i}\tilde{\chi_{i}}c\hat{a}_{i}}=$
$=\prod_{i}\left(1+\chi_{i}c\hat{a}_{i}^{\dagger}+\frac{(\chi_{i}c\hat{a}_{i}^{\dagger})^{2}}{2}\right)e^{\sum_{i}\xi_{0i}c\hat{a}_{i}^{\dagger}}\ket{0}\bra{0}e^{\sum_{i}\xi_{0i}c\hat{a}_{i}}\prod_{i}\left(1+\tilde{\chi_{i}}c\hat{a}_{i}+\frac{(\tilde{\chi_{i}}c\hat{a}_{i})^{2}}{2}\right).$

The creation operators $\hat{a}_{i}^{\dagger}$ that act on the input modes can
be written in terms of the operators $\hat{d}_{i}^{\dagger}$ that act on the
output modes:

$\hat{\mu}=\prod_{i}\left(1+\chi_{i}c\sum_{j}U_{ij}\hat{d}_{j}^{\dagger}+\frac{(\chi_{i}c\sum_{j}U_{ij}\hat{d}_{j}^{\dagger})^{2}}{2}\right)e^{\sum_{ij}\xi_{0i}cU_{ij}\hat{d}_{j}^{\dagger}}\ket{0}\cdot$
$\cdot\bra{0}e^{\sum_{ij}\xi_{0i}cU_{ij}^{*}\hat{d}_{j}}\prod_{i}\left(1+\tilde{\chi_{i}}c\sum_{j}U_{ij}^{*}\hat{d}_{j}+\frac{(\tilde{\chi_{i}}c\sum_{j}U_{ij}^{*}\hat{d}_{j})^{2}}{2}\right).$

We will denote

$\hat{\nu}(\vec{\xi}_{0}c)=e^{\sum_{ij}\xi_{0i}cU_{ij}\hat{d}_{j}^{\dagger}}\ket{0}\bra{0}e^{\sum_{ij}\xi_{0i}cU_{ij}^{*}\hat{d}_{j}}.$
(23)

We can expand the brackets in the expression for $\hat{\mu}$, leaving the
terms up to the second order:

$\prod_{i}\left(1+\chi_{i}c\sum_{j}U_{ij}\hat{d}^{\dagger}_{j}+\frac{(\chi_{i}c\sum_{j}U_{ij}\hat{d}^{\dagger}_{j})^{2}}{2}\right)=1+\sum_{j}\hat{d}^{\dagger}_{j}\sum_{i}c\chi_{i}U_{ij}+$
$+\sum_{jk}\hat{d}^{\dagger}_{j}\hat{d}^{\dagger}_{k}\left(\frac{1}{2}\sum_{i}c^{2}\chi^{2}_{i}U_{ij}U_{ik}+\sum_{i\neq
l}c\chi_{i}c_{l}\chi_{l}U_{ij}U_{lk}\right).$

$\prod_{i}\left(1+\tilde{\chi_{i}}c\sum_{j}U_{ij}^{*}\hat{d}_{j}+\frac{(\tilde{\chi_{i}}c\sum_{j}U_{ij}^{*}\hat{d}_{j})^{2}}{2}\right)=1+\sum_{j}\hat{d}_{j}\sum_{i}c\tilde{\chi_{i}}U_{ij}^{*}+$
$+\sum_{jk}\hat{d}_{j}\hat{d_{k}}\left(\frac{1}{2}\sum_{i}c^{2}\tilde{\chi}^{2}_{i}U_{ij}^{*}U_{ik}^{*}+\sum_{i\neq
l}c\tilde{\chi_{i}}c_{l}\tilde{\chi_{l}}U_{ij}^{*}U_{lk}^{*}\right).$

When we take the product of these two expressions, most of the resulting terms
will have zero expected value because of the properties of the normal
distribution. Then

$\hat{\mu}=\hat{\nu}(\vec{\xi}_{0}c)+\frac{1}{2}\sum_{i}\chi_{i}^{2}c^{2}\sum_{jk}U_{ij}U_{ik}\cdot\hat{d}^{\dagger}_{j}\hat{d}^{\dagger}_{k}\hat{\nu}(\vec{\xi}_{0}c)+\frac{1}{2}\sum_{i}\tilde{\chi_{i}}^{2}c^{2}\sum_{jk}U_{ij}^{*}U_{ik}^{*}\cdot\hat{\nu}(\vec{\xi}_{0}c)\hat{d}_{j}\hat{d}_{k}+$
$+\sum_{i}\chi_{i}\tilde{\chi_{i}}c^{2}\sum_{jk}U_{ij}U_{ik}^{*}\cdot\hat{d}^{\dagger}_{j}\hat{\nu}(\vec{\xi}_{0}c)\hat{d}_{k}+$
$+\frac{1}{4}\sum_{ij}\chi_{i}^{2}\tilde{\chi_{j}}^{2}c^{4}\sum_{klmn}U_{ik}U_{il}U_{jm}^{*}U_{jn}^{*}\cdot\hat{d}^{\dagger}_{k}\hat{d}^{\dagger}_{l}\hat{\nu}(\vec{\xi}_{0}c)\hat{d}_{m}\hat{d}_{n}+$
$+\sum_{i\neq
j}\chi_{i}\chi_{j}\tilde{\chi_{i}}\tilde{\chi_{j}}c^{4}\sum_{klmn}U_{ik}U_{jl}U_{im}^{*}U_{jn}^{*}\cdot\hat{d}^{\dagger}_{k}\hat{d}^{\dagger}_{l}\hat{\nu}(\vec{\xi}_{0}c)\hat{d}_{m}\hat{d}_{n}$

The integrals over $\chi_{i}$, $\tilde{\chi}_{i}$ result in specific moments
of the distribution, and the integral over $\xi_{0i}$ can be calculated using
Monte-Carlo methods. The final expression is

$Tr\left\\{\hat{\rho}_{out}\hat{\vec{n}}\right\\}=\frac{(\det\Sigma)^{N/2}}{\alpha^{N}(2\pi)^{3N/2}(\det\Gamma)^{N/2}}\cdot\int
d\vec{\xi}_{0}\Bigg{[}Tr\left\\{\hat{\nu}(\vec{\xi}_{0}c)\hat{\vec{n}}\right\\}+$
$+\frac{1}{2}\overline{\chi^{2}}c^{2}\sum_{ijk}U_{ij}U_{ik}\cdot
Tr\left\\{\hat{d}^{\dagger}_{j}\hat{d}^{\dagger}_{k}\hat{\nu}(\vec{\xi}_{0}c)\hat{\vec{n}}\right\\}+\frac{1}{2}\overline{\tilde{\chi}^{2}}c^{2}\sum_{ijk}U_{ij}^{*}U_{ik}^{*}\cdot
Tr\left\\{\hat{\nu}(\vec{\xi}_{0}c)\hat{d}_{j}\hat{d}_{k}\hat{\vec{n}}\right\\}+$
$+\overline{\chi\tilde{\chi}}c^{2}\sum_{ijk}U_{ij}U_{ik}^{*}\cdot
Tr\left\\{\hat{d}^{\dagger}_{j}\hat{\nu}(\vec{\xi}_{0}c)\hat{d}_{k}\hat{\vec{n}}\right\\}+$
$+\frac{1}{4}c^{4}\sum_{ij}\left(\left(\overline{\chi^{2}}\right)^{2}+2\delta_{ij}\left(\overline{\chi\tilde{\chi}}\right)^{2}\right)\sum_{klmn}U_{ik}U_{il}U_{jm}^{*}U_{jn}^{*}\cdot
Tr\left\\{\hat{d}^{\dagger}_{k}\hat{d}^{\dagger}_{l}\hat{\nu}(\vec{\xi}_{0}c)\hat{d}_{m}\hat{d}_{n}\hat{\vec{n}}\right\\}+$
$+\left(\overline{\chi\tilde{\chi}}\right)^{2}c^{4}\sum_{i\neq
j}\sum_{klmn}U_{ik}U_{jl}U_{im}^{*}U_{jn}^{*}\cdot
Tr\left\\{\hat{d}^{\dagger}_{k}\hat{d}^{\dagger}_{l}\hat{\nu}(\vec{\xi}_{0}c)\hat{d}_{m}\hat{d}_{n}\hat{\vec{n}}\right\\}\Bigg{]},$

where by Wick’s probability theorem
$\overline{\chi_{i}^{2}\tilde{\chi}_{j}^{2}}=\overline{\chi_{i}^{2}}\cdot\overline{\tilde{\chi}_{j}^{2}}+\overline{\chi_{i}\tilde{\chi}_{j}}\cdot\overline{\chi_{i}\tilde{\chi}_{j}}+\overline{\chi_{i}\tilde{\chi}_{j}}\cdot\overline{\chi_{i}\tilde{\chi}_{j}}=\left(\overline{\chi^{2}}\right)^{2}+2\delta_{ij}\left(\overline{\chi\tilde{\chi}}\right)^{2}$.

###  3.6 Calculating traces

In order to calculate $Tr\left\\{\hat{\rho}_{out}\hat{\vec{n}}\right\\}$, we
need to be able to calculate expression like
$Tr\left\\{\hat{\nu}(\vec{x})\hat{\vec{n}}\right\\}$,
$Tr\left\\{\hat{d}^{\dagger}_{j}\hat{d}^{\dagger}_{k}\hat{\nu}(\vec{x})\hat{\vec{n}}\right\\}$,
$Tr\left\\{\hat{\nu}(\vec{x})\hat{d}_{j}\hat{d}_{k}\hat{\vec{n}}\right\\}$,
$Tr\left\\{\hat{d}^{\dagger}_{j}\hat{\nu}(\vec{x})\hat{d}_{k}\hat{\vec{n}}\right\\}$,
etc, for different $\vec{x}$. The first one can be calculated fairly easily:

$Tr\left\\{\hat{\nu}(\vec{x})\hat{\vec{n}}\right\\}=Tr\left\\{e^{\sum_{ij}x_{i}U_{ij}\hat{d}_{j}^{\dagger}}\ket{0}\bra{0}e^{\sum_{ij}x_{i}U_{ij}^{*}\hat{d}_{j}}\ket{\vec{n}}\bra{\vec{n}}\right\\}=$
$=\bra{0}e^{\sum_{ij}x_{i}U_{ij}^{*}\hat{d}_{j}}\ket{\vec{n}}\bra{\vec{n}}e^{\sum_{ij}x_{i}cU_{ij}\hat{d}_{j}^{\dagger}}\ket{0}=$
$=\prod_{j}\bra{0}e^{\sum_{i}x_{i}U_{ij}^{*}\hat{d}_{j}}\ket{n_{j}}\bra{n_{j}}e^{\sum_{i}x_{i}U_{ij}\hat{d}_{j}^{\dagger}}\ket{0}=$
$=\prod_{j}\bra{0}\frac{\left(\sum_{i}x_{i}U_{ij}^{*}\hat{d}_{j}\right)^{n_{j}}}{n_{j}!}\ket{n_{j}}\bra{n_{j}}\frac{\left(\sum_{i}x_{i}U_{ij}\hat{d}_{j}^{\dagger}\right)^{n_{j}}}{n_{j}!}\ket{0}=$
$=\prod_{j}\left[\frac{\left(\sum_{i}x_{i}U_{ij}^{*}\right)^{n_{j}}}{\sqrt{n_{j}!}}\right]\cdot\left[\frac{\left(\sum_{i}x_{i}U_{ij}\right)^{n_{j}}}{\sqrt{n_{j}!}}\right]=$
$=\prod_{j}\frac{1}{n_{j}!}\left|\sum_{i}x_{i}U_{ij}\right|^{2n_{j}}.$

Now suppose we need to calculate
$Tr\left\\{(\hat{d}^{\dagger}_{1})^{q_{1}}...(\hat{d}^{\dagger}_{N})^{q_{N}}\hat{\nu}(\vec{x})(\hat{d}_{1})^{q_{1}}...(\hat{d}_{N})^{q_{N}}\hat{\vec{n}}\right\\}$.
First, we note that

$Tr\left\\{(\hat{d}^{\dagger}_{1})^{q_{1}}...(\hat{d}^{\dagger}_{N})^{q_{N}}\hat{\nu}(\vec{x})(\hat{d}_{1})^{p_{1}}...(\hat{d}_{N})^{p_{N}}\hat{\vec{n}}\right\\}=Tr\left\\{\hat{\nu}(\vec{x})(\hat{d}_{1})^{p_{1}}...(\hat{d}_{N})^{p_{N}}\hat{\vec{n}}(\hat{d}^{\dagger}_{1})^{q_{1}}...(\hat{d}^{\dagger}_{N})^{q_{N}}\right\\}=$
$=Tr\left\\{\hat{\nu}(\vec{x})\ket{\vec{n}-\vec{p}}\bra{\vec{n}-\vec{q}}\right\\}\cdot\prod_{j}\sqrt{n_{j}(n_{j}-1)...(n_{j}-p_{j}+1)n_{j}(n_{j}-1)...(n_{j}-q_{j}+1)}=$
$=Tr\left\\{\hat{\nu}(\vec{x})\ket{\vec{n}-\vec{p}}\bra{\vec{n}-\vec{q}}\right\\}\cdot\prod_{j}\sqrt{\frac{n_{j}!}{(n_{j}-p_{j})!}\frac{n_{j}!}{(n_{j}-q_{j})!}},$

where by e.g. $\ket{\vec{n}-\vec{p}}$ we mean
$\bigotimes_{i}\ket{n_{i}-p_{i}}$.

$=Tr\left\\{\hat{\nu}(\vec{x})\ket{\vec{n}-\vec{p}}\bra{\vec{n}-\vec{q}}\right\\}=\prod_{j}\bra{0}e^{\sum_{i}x_{i}U_{ij}^{*}\hat{d}_{j}}\ket{n_{j}-p_{j}}\bra{n_{j}-q_{j}}e^{\sum_{i}x_{i}U_{ij}\hat{d}_{j}^{\dagger}}\ket{0}=$
$=\prod_{j}\left[\frac{\left(\sum_{i}x_{i}U_{ij}^{*}\right)^{n_{j}-p_{j}}}{\sqrt{(n_{j}-p_{j})!}}\right]\cdot\left[\frac{\left(\sum_{i}x_{i}U_{ij}\right)^{n_{j}-q_{j}}}{\sqrt{(n_{j}-q_{j})!}}\right]=$
$=\prod_{j}\frac{1}{\sqrt{(n_{j}-p_{j})!(n_{j}-q_{j})!}}\cdot\frac{\left|\sum_{i}x_{i}U_{ij}\right|^{2n_{j}}}{\left(\sum_{i}x_{i}U_{ij}^{*}\right)^{p_{j}}\left(\sum_{i}x_{i}U_{ij}\right)^{q_{j}}}=$
$=\prod_{j}\frac{\left|\sum_{i}x_{i}U_{ij}\right|^{2n_{j}}}{\left(\sum_{i}x_{i}U_{ij}^{*}\right)^{p_{j}}\left(\sum_{i}x_{i}U_{ij}\right)^{q_{j}}}\cdot\frac{\sqrt{\frac{n_{j}!}{(n_{j}-p_{j})!}\frac{n_{j}!}{(n_{j}-q_{j})!}}}{n_{j}!}=$
$=Tr\left\\{\hat{\nu}(\vec{x})\hat{\vec{n}}\right\\}\prod_{j}\sqrt{\frac{n_{j}!}{(n_{j}-p_{j})!}\frac{n_{j}!}{(n_{j}-q_{j})!}}\frac{1}{\left(\sum_{i}x_{i}U_{ij}^{*}\right)^{p_{j}}\left(\sum_{i}x_{i}U_{ij}\right)^{q_{j}}}.$

Finally, we can write

$Tr\left\\{(\hat{d}^{\dagger}_{1})^{q_{1}}...(\hat{d}^{\dagger}_{N})^{q_{N}}\hat{\nu}(\vec{x})(\hat{d}_{1})^{p_{1}}...(\hat{d}_{N})^{p_{N}}\hat{\vec{n}}\right\\}=$
$=Tr\left\\{\hat{\nu}(\vec{x})\hat{\vec{n}}\right\\}\prod_{j}\frac{n_{j}!}{(n_{j}-p_{j})!}\frac{n_{j}!}{(n_{j}-q_{j})!}\frac{1}{\left(\sum_{i}x_{i}U_{ij}^{*}\right)^{p_{j}}\left(\sum_{i}x_{i}U_{ij}\right)^{q_{j}}}.$

##  4 Algorithm overview

The goal of the algorithm is to calculate the probability of a state
$\ket{\vec{n}}$, given $\vec{n}$, $\alpha$, $c$, $s$ and $U$. We assume that
the Taylor series expansion is done up to the desired order before computation
starts. The integrals over $\chi_{i}$ and $\tilde{\chi}_{i}$ should also be
computed (it can be done analytically via Wick’s probability theorem).

We start by calculating two-variable covariance matrix $\Sigma$ using $\alpha$
and $s$. We now select $\Gamma$ in the way specified above such that it
minimizes the series expansion parameter $\varepsilon$. In order to compute
the integrals over $\xi_{0i}$, we sample $\xi_{0i}$ for each $i$ from a normal
distribution $\mathbb{N}(0,\overline{\xi_{0}^{2}})$.

We now compute $Tr\left\\{\hat{\mu}\hat{\vec{n}}\right\\}$, which by linearity
consists in computing traces of the form described above; for each sample
$\vec{\xi}_{0}$ we need only a polynomial number of operations.

Finally, we take an average over our samples and multiply by the necessary
constant terms.

##  5 Taylor series convergence for actual experimental conditions

We have discussed above the fact that the role of the ”perturbation parameter”
in the series expansion is played by
$c^{2}\cdot\min(\overline{\chi^{2}},|h|)$, which we can choose to be equal to
$\varepsilon=\frac{1}{2}\frac{c^{2}}{1/\alpha+s^{2}}$. This parameter depends
on the experimental conditions (i.e. the squeezing parameter of the input
state $\alpha$ and loss level $s^{2}$). The smaller this parameter is, the
faster the series will converge. Thus, the best conditions for this algorithm
are achieved when the loss level $s^{2}$ is high and the squeezing parameter
$\alpha$ is low. Let us consider actual experimental implementation of the
gaussian boson sampling problem, and estimate how small this parameter is in
those conditions.

Let’s consider the relation between $\alpha$ and the average amount of photons
per state $\langle n\rangle$. If the squeezing parameter is
$\zeta=re^{i\varphi}$, then $\alpha=\tanh{r}$, while $\langle
n\rangle=\sinh^{2}{r}$.

In a paper by Zhong et al. [8] 25 PPKTP crystals were used to produce 25 two-
mode squeezed states, which is equivalent to 50 single-mode squeezed states.
The average number of photons registered by the detectors is 43. Thus, the
average amount of photons per mode $\langle n\rangle$ is around
$\frac{43}{50}$; $r=\operatorname{arcsinh}(\sqrt{\langle n\rangle})\approx
0.855$, $\alpha=\tanh{r}\approx 0.694$. The average collection efficiency is
said to be $c^{2}=0.628$. Then,
$\varepsilon=\frac{1}{2}\frac{c^{2}}{\frac{1}{\alpha}+s^{2}}\approx 0.18$.

In another paper by Zhong et al. [9], the average amount of photons produced
was increased to $70$ at maximum pump intensity. This corresponds to
$\alpha\approx 0.76$. The overall transmission rate in the experiment is said
in the paper to be $48\%$ and $54\%$ for different settings, so we take
$s^{2}\approx 0.5$. This yields $\varepsilon\approx 0.14$.

The conclusion that we draw is that even in large GBS experiments which are
said to demonstrate quantum advantage the conditions are such that
$\varepsilon$ is fairly small, and not many orders of the series need to be
calculated to produce an approximation.

##  6 Implementation details

###  6.1 Contraction precomputation

Let’s consider the term

$\frac{1}{2}\overline{\chi^{2}}c^{2}\sum_{ijk}U_{ij}U_{ik}\cdot
Tr\left\\{\hat{d}^{\dagger}_{j}\hat{d}^{\dagger}_{k}\hat{\nu}(\vec{\xi}_{0}c)\hat{\vec{n}}\right\\}.$

We can rewrite it as

$\frac{1}{2}\overline{\chi^{2}}c^{2}\sum_{jk}Tr\left\\{\hat{d}^{\dagger}_{j}\hat{d}^{\dagger}_{k}\hat{\nu}(\vec{\xi}_{0}c)\hat{\vec{n}}\right\\}\sum_{i}U_{ij}U_{ik}=\frac{1}{2}\overline{\chi^{2}}c^{2}\sum_{jk}Tr\left\\{\hat{d}^{\dagger}_{j}\hat{d}^{\dagger}_{k}\hat{\nu}(\vec{\xi}_{0}c)\hat{\vec{n}}\right\\}T_{jk},$

where $T_{jk}=\sum_{i}U_{ij}U_{ik}$ is a contraction of $U$ with itself. It
depends only on $U$, and can be calculated before sampling $\vec{\xi}_{0}$,
which reduces the amount of operations required to calculate each probability
sample from a $\vec{\xi}_{0}$ sample.

###  6.2 Factorial fractions precomputation

In calculating traces of the form described above, we need to calculate
factorial fractions of the form $\frac{m!}{(m-p)!}\equiv F^{m}_{p}$, where
$0\leq p\leq m$. Since the target state $\hat{\vec{n}}$ is fixed,
$m\leq\max(n_{i})$.

###  6.3 Reusing $\sum_{i}x_{i}U_{ij}$

During calculation, while calculating each trace, we can calculate
$\sum_{i}x_{i}U_{ij}$ only once for each $\vec{\xi}_{0}$ sample and then reuse
it, thus using less operations to calculate each trace. Let’s denote
$S_{j}=\sum_{i}x_{i}U_{ij}$; $\vec{S}=U^{T}\vec{x}$. Then,

$Tr\left\\{\hat{\nu}(\vec{x})\hat{\vec{n}}\right\\}=\prod_{j}\frac{1}{n_{j}!}\left|S_{j}\right|^{2n_{j}}$
(24)

and

$Tr\left\\{(\hat{d}^{\dagger}_{1})^{q_{1}}...(\hat{d}^{\dagger}_{N})^{q_{N}}\hat{\nu}(\vec{x})(\hat{d}_{1})^{p_{1}}...(\hat{d}_{N})^{p_{N}}\hat{\vec{n}}\right\\}=$
$=Tr\left\\{\hat{\nu}(\vec{x})\hat{\vec{n}}\right\\}\prod_{j}\frac{n_{j}(n_{j}-1)...(n_{j}-p_{j}+1)n_{j}(n_{j}-1)...(n_{j}-q_{j}+1)}{\left(S_{j}^{*}\right)^{p_{j}}\left(S_{j}\right)^{q_{j}}}.$

##  7 Complexity analysis

###  7.1 Precomputation

In this section we will analyze the computational complexity of
precomputation. By precomputation we mean the calculations that need to be
carried out only once before $\vec{\xi}_{0}$ sampling and before calculating
probability samples for each $\vec{\xi}_{0}$. The multiplicative constant
before the integral sign
$\frac{(\det\Sigma)^{N/2}}{\alpha^{N}(2\pi)^{3N/2}(\det\Gamma)^{N/2}}$ can be
calculated with $O(N)$ multiplication operations. For each term in the
resulting sum, we will define its order to be the number of variables $\chi$
and $\tilde{\chi}$, or, equivalently, the power of the loss parameter $c$.
Thus, the term

$\overline{\chi\tilde{\chi}}c^{2}\sum_{ijk}U_{ij}U_{ik}^{*}\cdot
Tr\left\\{\hat{d}^{\dagger}_{j}\hat{\nu}(\vec{\xi}_{0}c)\hat{d}_{k}\hat{\vec{n}}\right\\}$

will be of the second order. Then, each term of the order $K$ will have a
contraction of the form

$\sum_{j_{1}...j_{K}}U_{i_{1}j_{1}}U_{i_{2}j_{2}}...U_{i_{K}j_{K}}$ (25)

where some of the $U_{ji}$ can be conjugated. This leaves at most $K+1$
different ways to conjugate the factors. Each contraction has $K$ free
indices, and calculating the sum requires $N^{K}$ additions and $N^{K}(K-1)$
multiplications. The total number of additions is $N^{2K}$ and the number of
multiplications is $N^{2K}(K-1)$, where $K$ is the maximum order we choose to
calculate.

Calculating all $F^{m}_{p}\equiv\frac{m!}{(m-p)!}$ for $0\leq p\leq
m\leq\max(n_{i})$ requires only around $\frac{\max(n_{i})^{2}}{2}$
multiplications, since $\forall m~{}~{}F^{m}_{0}=1$, $F^{m}_{1}=m$,
$F^{m}_{2}=m(m-1)=(m-1)F^{m}_{1}$, …, $F^{m}_{k}=(m-k+1)F^{m}_{k-1}$.

###  7.2 Probability sample computation

Here we will analyze the computational complexity of calculating a single
probability sample given $\vec{\xi}_{0}$. We will assume that the terms are
calculated up to some order $K$.

Calculating the trace

$Tr\left\\{\hat{\nu}(\vec{x})\hat{\vec{n}}\right\\}=\prod_{j}\frac{1}{n_{j}!}\left|\sum_{i}x_{i}U_{ij}\right|^{2n_{j}}.$

requires one multiplication of an $N\times N$ matrix by a $N$-dimensional
vector, $N$ exponentiation operations and $2N$ multiplication operations. This
calculations needs to be done only once for each $\vec{x}$. Calculating any
other trace of the form

$Tr\left\\{(\hat{d}^{\dagger}_{1})^{q_{1}}...(\hat{d}^{\dagger}_{N})^{q_{N}}\hat{\nu}(\vec{x})(\hat{d}_{1})^{p_{1}}...(\hat{d}_{N})^{p_{N}}\hat{\vec{n}}\right\\}$

requires $2N$ exponentiation operations and $4N$ multiplication operations
(since factorial fractions are precomputed).

The number of terms for a given order $K$ is $N^{K}$ times the number of
different non-zero $K$-th order moments
$\overline{\chi_{i_{1}}\dotsm\chi_{i_{r}}\tilde{\chi}_{i_{r+1}}\dotsm\tilde{\chi}_{i_{K}}}$.
The exact amount is hard to calculate, but the total number of moments
(including those that are zero) is $(K+1)N^{K}$. Thus, the maximum amount of
terms required to is $(K+1)N^{2K}$.

Since the amount of operations required to calculate each term is $O(N)$, the
total computational complexity of calculating a probability sample for a given
$\vec{\xi}_{0}$ is $O\left(K\cdot N^{2K}\right)$

##  8 Results

Below are the results of probability calculation for $N=5$ for different
output states. The calculated probabilities are compared to exact solutions.
The parameters are: $\alpha=0.9$, $c=s=\sqrt{2}$. Number of samples is $4096$.

Figure 1: Probability calculation for $5$ modes for different $2$-photon
output states. Figure 2: Graph of the average probability and the standard
deviation calculated up to different orders for different numbers of samples.
The state for this graph is $2$-photon.

These results show that for calculating a single output state probability
accurately the number of samples needs to be on the order of $10^{4}$. Below
are the results of using fewer samples per state, but instead of comparing
individual probabilities we look at cosine similarity between exact and
approximated probability distributions over all 2-photon states.

Figure 3: Convergence of the cosine similarity between estimated probability
distribution over the set of all 2-photon states and ground truth for
different $N$.

The above graph suggests that the number of samples per state needed to
approximate the distribution doesn’t depend much on $N$. It is computationally
hard to check this when comparing to the exact solution, but if we assume that
the cosine similarity converges to a value close to $1$, we can estimate how
quickly it converges. Below we look at the cosine similarity between a
distribution calculated with $K$ samples per state and a distribution
calculated with $K+10$ samples per state for different K. Figure 4 suggests
more strongly that the number of samples per state required for accurate
approximation is not really influenced by $N$

Figure 4: Cosine similarity between probability distribution over the set of
all 2-photon states after $K$ samples and after $K+10$ samples for different
$N$.

Below are benchmark results that show average precomputation time, which
depends only on $N$, and time per sample, which depends on $N$ and the amount
of photons $M$ in the target state.

Figure 5: Precomputation time on an Intel i5 CPU in ms versus the number of
modes. Figure 6: Average time per sample on an Intel i5 CPU versus the number
of modes for states with different photon numbers.

These results show that even $N=40$ mode GBS can be simulated on an average
laptop using this algorithm.

##  9 Conclusion

In this paper we have presented a new algorithm for the approximate
calculation of the probability of observing a given output state in Gaussian
boson sampling instance. We have discussed various implementation details that
help to reduce the number of operations needed to calculate each probability
sample. We also analyze the total computational complexity both of the
calculations that need to be carried out once for each specific problem and of
computing each probability sample.

This algorithm relies on the Taylor series expansion where the ”perturbation”
parameter is dependent on the problem conditions. The algorithm consists in
calculating the terms of this Taylor series up to some finite order. For a
fixed maximum order, the computational complexity of the algorithm is
polynomial in $N$.

We have demonstrated that increasing the maximum order does increase the
accuracy of the answer. We have also measured precomputation and sampling time
for a regular CPU, showing that even large instances of Gaussian boson
sampling ($N\approx 40$) can be solved in reasonable time.

We have considered recent GBS experiments and estimated the parameters of the
problem for those conditions. We believe that in those conditions our
algorithm doesn’t require many orders of the Taylor series to be calculated
for approximating a probability of an output state.

## References

  * [1] Peter W Shor  “Polynomial-Time Algorithms for Prime Factorization and Discrete Logarithms on a Quantum Computer”  In _SIAM Journal on Computing_ 26.5 Society for Industrial & Applied Mathematics (SIAM), 1997, pp. 1484–1509  DOI: 10.1137/s0097539795293172
  * [2] A. P. Lund et al.  “Boson Sampling from a Gaussian State”  In _Physical Review Letters_ 113.10 American Physical Society (APS), 2014  DOI: 10.1103/physrevlett.113.100502
  * [3] Bryan T. Gard et al.  “An Introduction to Boson-Sampling”  In _From Atomic to Mesoscale_ WORLD SCIENTIFIC, 2015, pp. 167–192  DOI: 10.1142/9789814678704˙0008
  * [4] Marco Bentivegna et al.  “Experimental scattershot boson sampling”  In _Science Advances_ 1.3, 2015, pp. e1400255  DOI: 10.1126/sciadv.1400255
  * [5] Craig S. Hamilton et al.  “Gaussian Boson Sampling”  In _Physical Review Letters_ 119.17 American Physical Society (APS), 2017  DOI: 10.1103/physrevlett.119.170501
  * [6] A.. Lund et al.  “Boson Sampling from a Gaussian State”  In _Phys. Rev. Lett._ 113 American Physical Society, 2014, pp. 100502  DOI: 10.1103/PhysRevLett.113.100502
  * [7] Han-Sen Zhong et al.  “Experimental Gaussian Boson sampling”  In _Science Bulletin_ 64.8 Elsevier BV, 2019, pp. 511–515  DOI: 10.1016/j.scib.2019.04.007
  * [8] Han-Sen Zhong et al.  “Quantum computational advantage using photons”  In _Science_ 370.6523 American Association for the Advancement of Science (AAAS), 2020, pp. 1460–1463  DOI: 10.1126/science.abe8770
  * [9] Han-Sen Zhong et al.  “Phase-Programmable Gaussian Boson Sampling Using Stimulated Squeezed Light”  In _Physical Review Letters_ 127.18 American Physical Society (APS), 2021  DOI: 10.1103/physrevlett.127.180502
  * [10] Changhun Oh, Kyungjoo Noh, Bill Fefferman and Liang Jiang  “Classical simulation of lossy boson sampling using matrix product operators”  In _Physical Review A_ 104.2 American Physical Society (APS), 2021  DOI: 10.1103/physreva.104.022407
  * [11] Raúl García-Patrón, Jelmer J. Renema and Valery Shchesnovich  “Simulating boson sampling in lossy architectures”  In _Quantum_ 3 Verein zur Forderung des Open Access Publizierens in den Quantenwissenschaften, 2019, pp. 169  DOI: 10.22331/q-2019-08-05-169
  * [12] A.. Popova and A.. Rubtsov  “Cracking the Quantum Advantage threshold for Gaussian Boson Sampling”  arXiv, 2021  DOI: 10.48550/ARXIV.2106.01445
  * [13] Haoyu Qi, Daniel J. Brod, Nicolá s Quesada and Raúl García-Patrón  “Regimes of Classical Simulability for Noisy Gaussian Boson Sampling”  In _Physical Review Letters_ 124.10 American Physical Society (APS), 2020  DOI: 10.1103/physrevlett.124.100502
  * [14] Scott Aaronson and Daniel J. Brod  “BosonSampling with lost photons”  In _Physical Review A_ 93.1 American Physical Society (APS), 2016  DOI: 10.1103/physreva.93.012335
  * [15] Changhun Oh et al.  “Classical algorithm for simulating experimental Gaussian boson sampling”, 2023  arXiv:2306.03709 [quant-ph]
  * [16] J. Hubbard  “Calculation of Partition Functions”  In _Phys. Rev. Lett._ 3 American Physical Society, 1959, pp. 77–78  DOI: 10.1103/PhysRevLett.3.77
  * [17] R.. Stratonovich  “On a Method of Calculating Quantum Distribution Functions”  In _Soviet Physics Doklady_ 2, 1957, pp. 416


# Crystal-field effects in the formation of Wigner-molecule supercrystals in
moiré TMD superlattices

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

(27 May 2024)

###### Abstract

For moiré bilayer TMD superlattices, full-configuration-interaction (FCI)
calculations are presented that take into account both the intra-moiré-
quantum-dot (MQD) charge-carrier Coulombic interactions, as well as the
crystal-field effect from the surrounding moiré pockets (inter-moiré-QD
interactions). Such FCI calculations enable an effective computational
embedding strategy and allow for a complete interpretation of the
counterintuitive experimental observations reported recently in the context of
moiré TMD superlattices at integer fillings $\nu=2$ and 4. Two novel states of
matter are reported: (i) a genuinely quantum-mechanical supercrystal of
sliding Wigner molecules (WMs) for unstrained moiré TMD materials (when the
crystal field is commensurate with the trilobal symmetry of the confining
potential in each embedded MQD) and (ii) a supercrystal of pinned Wigner
molecules when strain is involved and the crystal field is incommensurate with
the trilobal symmetry of the confining potential in each embedded MQD. The
case of $\nu=3$ is an exception, in that both unstrained and strained cases
produce a supercrystal of pinned WMs, which is due to the congruence of
intrinsic (that of the WM) and external (that of the confining potential of
the MQD) $C_{3}$ point-group symmetries. Furthermore, it is shown that the
unrestricted Hartree-Fock approach fails to describe the supercrystal of
sliding WMs in the unstrained case, providing a qualitative agreement only in
the case of a supercrystal of pinned WMs.

##  I Introduction

Earlier studies have revealed a novel fundamental-physics aspect in artificial
twodimensional (2D) nanosystems in the regime of strong interelectronic
correlations, namely formation of quantum Wigner molecules (WMs), originally
described theoretically [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,
16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26] in 2D semiconductor (parabolic,
elliptic, and double-well) quantum dors (QDs), as well as in trapped ultracold
atoms, and subsequently observed experimentally in GaAs QDs [27, 28, 29, 30],
Si/SiGe QDs [31], planar germanium QDs [32] and 1D carbon-nanotubes [33].

Remarkably, most recent work [34, 35, 36] has extended the WM portfolio to the
newly emerging and highly regarded (due to the potential for fundamental-
physics discoveries and for advancing quantum-device applications) field of
transition-metal dichalcogenide (TMD) moiré materials and superlattices.
Specifically, TMD moiré superlattices are highly pursued for solving the
question of scalability in quantum computer architectures.

Figure 1: FCI ground-state (with total spin $S=0$) charge densities for $N=2$
holes in the embedded central MQD (red triangle in the schematics) taking into
consideration the crystal-field effect generated by the charge carriers in the
surrounding six moiré pockets (see blue hexagon, with radius $a_{M}$, in the
schematics). (a) and (c) $\kappa=5$ in agreement with the case of a TMD moiré
superlattice in hBN environment. (b) and (d) were calculated, for purpose of
comparison, with $\kappa=2$, corresponding to a stronger Coulomb repulsion. In
(a) and (b), all six surrounding MQDs (see schematics) were populated with
$Q=Ne$ charge carriers. In (c) and (d), only one of the surrounding MQDs (the
one to the right, see schematics) has been populated with $Q=Ne$ charge
carriers. The charge-carrier distributions in each surrounding MQD have been
represented by three point charges $Q/3$ placed at the apices of an
equilateral triangle (see schematics), with radius $a_{M}/6$. Remaining
parameters: effective mass $m^{*}=0.90m_{e}$, moiré lattice constant
$a_{M}=9.8$ nm, depth of a moiré pocket $v_{0}=10.3$ eV, and
$\phi=20^{\circ}$. Charge densities are in units of 1/nm2.

Building on the demonstrated formation [34] of WMs in the isolated moiré
pockets [most often referred to as moiré quantum dots (MQDs)] at integer
fillings, we investigate here the effects arising from embedding such single
MQDs in the moiré superlattice structure, which is the actual system that was
addressed experimentally [36]. Specifically, this paper adresses the key
question whether the experimental observations in the unstrained-lattice case
do reflect formation of a collective WM supercrystal, i.e., a novel, more
complex insulating phase of matter beyond the well-known generalized Wigner
crystal, or the trivial case of a collection of independent MQDs. The present
investigations (see below) support formation of a WM supercrystal 111We note
that the experimental results for the unstrained case in Ref. [36] were
necessarily interpreted with formation of sliding WMs in isolated MQDs (for
$N=2$ confined charge carriers, see Refs. [34, 36]; for $N=4$ confined charge
carriers, see Ref. [34]), given the fact that the mean-field Hartree-Fock
approach (that allows calculations for a large number of charge carriers
distributed over several moiré pockets) is unreliable [35, 36] for describing
sliding WMs. having sliding WMs as building blocks for $\nu=2$ or $\nu=4$,
with the filling $\nu=3$ being an exception having a pinned WM as a building
block.

Specifically, this paper focuses on the effects on WM formation in a given MQD
resulting from the crystal-field generated by the charge carriers (of same
nature) in the surrounding MQDs (moiré pockets) located at the apices of a
regular hexagon in the TMD moiré superlattice. This approach is inspired by
the well-known crystal-field theory [38, 39] in molecular physics, and is
applicable here because the moiré lattice constant $a_{M}$ is much larger than
the extent of the Wigner molecule within a single MQD. Two different
computational methodologies will be used in this endeavor, namely: (i) the
full configuration interaction (FCI) [40, 8, 9, 41, 42, 14, 17, 43, 21, 22,
26, 26, 44, 45] and (ii) the spin-and-space unrestricted Hartree-Fock (sS-UHF)
[1, 46, 5, 14, 47].

##  II Many-body Hamiltonian (including crystal field from surrounding moiré
pockets)

Following earlier established literature [48, 49, 50, 34, 36], we approximate
the potential of the 2D TMD moiré superlattice that confines the extra charge
carriers as

$\displaystyle V({\bf r})=-2v_{0}\sum_{i=1}^{3}\cos({\bf G}_{i}\cdot{\bf
r}+\phi),$ (1)

where ${\bf G}_{i}=[(4\pi/\sqrt{3}a_{M})(\sin(2\pi i/3),\cos(2\pi i/3))]$ are
the moiré reciprocal lattice vectors. The materials specific parameters of
$V({\bf r})$ are $v_{0}$ (which can also be experimentally controlled through
voltage biasing), the moiré lattice constant $a_{M}$, and the angle $\phi$.
$a_{M}$ is typically of the order of 10 nm, which is much larger than the
lattice constant of the monolayer TMD material (typically a few Å).

The parameter $\phi$ controls the strength of the trilobal $C_{3}$ anisotropy
in each MQD potential pocket. This trilobal anisotropy can be seen by
expanding $V({\bf r})$ in Eq. (1) in powers of $r$, and defining an
approximate confining potential, $V_{\rm MQD}({\bf r})$, for a single MQD as
follows:

$\displaystyle\begin{split}V_{\rm MQD}({\bf r})\equiv V({\bf
r})+6v_{0}\cos(\phi)\approx m^{*}\omega_{0}^{2}r^{2}/2+{\cal
C}\sin(3\theta)r^{3}.\end{split}$ (2)

with $m^{*}\omega_{0}^{2}=16\pi^{2}v_{0}\cos(\phi)/a_{M}^{2}$, and ${\cal
C}=16\pi^{3}v_{0}\sin(\phi)/(3\sqrt{3}a_{M}^{3})$; $m^{*}$ is the effective
mass and the expansion of $V({\bf r})$ can be restricted to the terms up to
$r^{3}$. $(r,\theta)$ are the polar coordinates of the position vector ${\bf
r}$.

The effective many-body Hamiltonian associated with a given MQD embedded in
the moiré superlattice is given by

$\displaystyle H_{\rm MB}=H_{\rm MQD}+H_{\rm CF},$ (3)

where

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

with $\kappa$ being the dielectric constant, and

$\displaystyle H_{\rm CF}=\sum_{i=1}^{N}\sum_{j=1}^{3{\cal
M}}\frac{eQ}{3\kappa|{\bf r}_{i}-{\bf R}_{j}|}$ (5)

is the crystal-field potential from the surrounding moiré pockets, ${\cal M}$
denoting the number of surrounding moiré pockets populated with charge
carriers, and the ${\bf R}_{j}$’s being the positions of the point charges
that mimick the charge carrier distributions. In each one of the surrounding
moiré pockets, we use three point charges at the apices of an equilateral
triangle with a total charge $Q=Ne$, thus respecting the $C_{3}$ trilobal
symmetry; see solid dots in the schematics in the figures.

A brief outline of the FCI and sS-UHF methodologies, used to solve the
corresponding many-body Schrödinger equation, is presented in the Appendices.
Before proceeding with the computational results, we point out that the
arrangement having the equilateral three-point-charge distributions in all the
six surrounding moiré pockets [see schematic in Fig. 1(a)] displays an overall
$C_{6}$ point-group symmetry, which is commensurate with the trilobal $C_{3}$
symmetry associated with the potential confinement [Eq. (2)] of a single MQD.
Below, we will investigate the effect of both such a commensurate crystal
field, as well as the effect of an incommensurate crystal field (IcCF)
generated by placing the three-point-charge distribution in only one of the
six surrounding moiré pockets [see schematic in Fig. 1(c)].

Figure 2: sS-UHF ground-state (with total-spin projection $S_{z}=0$) charge
densities for $N=2$ holes in the central MQD (red triangle in the schematics)
taking into consideration the crystal-field effect generated by the charge
carriers in the surrounding six moiré pockets (see blue hexagon, with radius
$a_{M}$, in the schematics). (a) and (c) $\kappa=5$ in agreement with the case
of a TMD moiré superlattice in hBN environment. (b) and (d) were calculated,
for purpose of comparison, with $\kappa=2$, corresponding to a stronger
Coulomb repulsion. In (a) and (b), all six surrounding MQDs (see schematics)
were populated with $Q=Ne$ charge carriers. In (c) and (d), only one of the
surrounding MQDs (the one to the right, see schematics) has been populated
with $Q=Ne$ charge carriers. The charge-carrier distributions in each
surrounding MQD have been represented by three point charges $Q/3$ placed at
the apices of an equilateral triangle (see schematics), with radius $a_{M}/6$.
Remaining parameters: effective mass $m^{*}=0.90m_{e}$, moiré lattice constant
$a_{M}=9.8$ nm, depth of a moiré pocket $v_{0}=10.3$ eV, and
$\phi=20^{\circ}$. Charge densities are in units of 1/nm2.

##  III Results for filling factor $\nu=2$

The FCI and sS-UHF charge densities (CDs) for a moiré TMD (e.g., WS2 [36])
superlattice at integer filling factor $\nu=2$ calculated using the many-body
Hamiltonian (…) are displayed in Fig. 1 and Fig. 2, respectively. The
configuration of the equilateral three-point-charge distributions that mimick
the crystal-field effect from the surrounding six moiré pockets is illustrated
in the schematics drawn in the upper right corner of each frame; see figure
captions for details and parameters. The total charge associated with each
equilateral three-point-charge distribution is $Q=2e$. (Note that $Q/e$ can be
taken different than $\nu$ in the case of a noninteger filling.)

The FCI calculations [see Figs. 1(a) and 1(b)] demonstrate that the full
crystal-field effect (all six surrounding moiré pockets are contributing)
maintains the charge densities that are associated with sliding WMs. Namely,
the CDs preserve the $C_{3}$ symmetry of the confining single-moiré-pocket
potential, as was the case with the isolated MQDs [34]. Note that even a
stronger Coulomb repulsion ($\kappa=2$ instead of $\kappa=5$, which
corresponds to the experimental setup) does not destroy the the $C_{3}$
symmetry; see Fig. 1(b). This symmetry-preserving behavior is counterintuitive
for a crystalline material in that the CDs necessarily exhibit three modest
humps instead of the naive expectation of a symmetry-breaking double hump
associated with two Wigner-crystal-type [51, 52, 53, 54] localized charge
carriers (one hump per charge carrier). Note further that this
counterintuitive behavior is in remarkable agreement with the recent
experimental findings [36] for the unstrained case at filling $\nu=2$.

Figure 3: Charge densities for $N=4$ holes in the central MQD (red triangle in
the schematics) taking into consideration the crystal-field effect generated
by the charge carriers in the surrounding six moiré pockets (see blue hexagon,
with radius $a_{M}$, in the schematics). (a) and (b) FCI ground-states (with
total spin $S=1$ and total-spin projection $S_{z}=0$). (c) and (d) sS-UHF
ground-states (with total spin projection $S_{z}=0$). Dielectric constant
$\kappa=5$ in agreement with the case of a TMD moiré superlattice in an hBN
environment. In (a) and (c), all six surrounding MQDs (see schematics) were
populated with $Q=Ne$ charge carriers. In (b) and (d), only one of the
surrounding MQDs (the one to the right, see schematics) has been populated
with $Q=Ne$ charge carriers. The charge-carrier distributions in each
surrounding MQD have been represented by three point charges $Q/3$ placed at
the apices of an equilateral triangle (see schematics), with radius $a_{M}/6$.
Remaining parameters: effective mass $m^{*}=0.90m_{e}$, moiré lattice constant
$a_{M}=9.8$ nm, depth of a moiré pocket $v_{0}=10.3$ eV, and
$\phi=20^{\circ}$. Charge densities are in units of 1/nm2.  Figure 4: Charge
densities for $N=3$ holes in the central MQD (red triangle in the schematics)
taking into consideration the crystal-field effect generated by the charge
carriers in the surrounding six moiré pockets (see blue hexagon, with radius
$a_{M}$, in the schematics). (a) FCI ground-states (with total spin $S=1/2$
and total-spin projection $S_{z}=1/2$). (b) sS-UHF ground-states (with total
spin projection $S_{z}=1/2$). Dielectric constant $\kappa=5$ in agreement with
the case of a TMD moiré superlattice in an hBN environment. In (a) and (b),
all six surrounding MQDs (see schematics) were populated with $Q=Ne$ charge
carriers. The charge-carrier distributions in each surrounding MQD have been
represented by three point charges $Q/3$ placed at the apices of an
equilateral triangle (see schematics), with radius $a_{M}/6$. Remaining
parameters: effective mass $m^{*}=0.90m_{e}$, moiré lattice constant
$a_{M}=9.8$ nm, depth of a moiré pocket $v_{0}=10.3$ eV, and
$\phi=20^{\circ}$. Charge densities are in units of 1/nm2.

The CDs in Figs. 2(a) and 2(b) demonstrate that the sS-UHF is unable to
properly describe the unstrained TMD moiré supercrystal at filling $\nu=2$.
Indeed, these sS-UHF CDs exhibit a pinned WM with two well-localized humps in
accordance with the traditional Wigner-crystal-type [51, 52, 53, 54]
expectation (one hump per charge carrier). We stress again that these sS-UHF
CDs do break the $C_{3}$ symmetry of the single-moiré-pocket confinement and
disagree with the experimental observations [36].

Heretofore, we used a crystal field with $C_{6}$ symmetry (associated with
that of the six surrounding moiré pockets), which is commensurate with the
$C_{3}$ symmetry of the confinement of the embedded MQD. A natural question
arising at this point concerns the effect of a crystal field that has a
symmetry non-commensurate with the $C_{3}$ symmetry of the single MQD.
Experimentally, such a non-commensurability was generated [36] by deforming
the moiré lattice through straining. Here, we generate a non-commensurability
situation by maintaining the equilateral three-point-charge configuration only
in a single moiré pocket [the one on the right; see schematics in Figs. 1(c,d)
and 2(c,d)].

From the FCI CDs in Figs. 1(c,d), it is apparent that the IcCF produces an
azimuthally pinned WM exhibiting two antipodal humps. These pinned-WM CDs do
break the $C_{3}$ symmetry of the single MQD. In addition, the extent of
localization for each charge carrier becomes more pronounced for stronger
Coulomb repulsion; compare the case of $\kappa=5$ [Fig. 1(c)] to $\kappa=2$
[Fig. 1(d)]. Concerning the corresponding sS-UHF CDs, it is obvious that, in
qualitative agreement with the FCI result, the IcCF CDs result also in pinned
two-humped WMs, exhibiting, nevertheless, a stronger charge-carrier
localization; contrast Fig. 1(c) to Fig. 2(c) and Fig. 1(d) to Fig. 2(d). We
mention here that our IcCF FCI CDs are in agreement with the recent
experimental observations [36] of two-humped pinned WMs per moiré pocket in
the case of strained moiré superlattices at a filling $\nu=2$.

##  IV Results for filling factor $\nu=4$

The ability of the crystal-field with the commensurate $C_{6}$ symmetry to
preserve, and even enhance, the formation of a quantum sliding WM in the
central MQD is revealed in an even more spectacular way in the case of the
filling factor $\nu=4$, where Ref. [36] has reported within each MQD the
observation of CDs with 3 humps ($C_{3}$ symmetry) for the unstrained lattice,
but with 4 humps (broken symmetry) for the strained lattice. Our FCI CDs that
take into account the crystal-field effect are in agreement with the
experiment; see the 3-hump, sliding-WM CD in Fig. 3(a) (commensurate crystal
field) and the 4-hump, distorted-pinne-WM CD in Fig. 3(b) (IcCF). We further
note that the sS-UHF fails to describe the sliding WM, yielding instead a
broken-$C_{3}$-symmetry CD (pinned WM) with 4 well defined humps even for the
case of the commensurate crystal field; see Fig. 3(c). In the case of an IcCF,
the sS-UHF agrees qualitatively with the FCI CD as it describes a deformed,
4-hump, pinned WM; compare Fig. 3(b) and Fig. 3(d), the charge-carrier
localization being more pronounced in the sS-UHF CD.

##  V Results for filling factor $\nu=3$

Compared to the cases of $\nu=2$ and $\nu=4$, the case of $\nu=3$ is
exceptional, a fact that was reported already in Ref. [34] where formation of
WMs was investigated for the case of an isolated MQD in the absence of a
crystal field from the surrounding moiré pockets. Specifically, the 3-hump
pinned WM in the central MQD has a $C_{3}$ symmetry which coincides with that
of the confining potential. This congruence of intrinsic and external point-
group symmetries prohibits the formation of a sliding WM and yields a strongly
pinned WM. The commensurate crystal field does not alter this CD distribution,
a fact that is reflected in the FCI calculation [see Fig. 4(a)]. In the
$\nu=3$ case, the sS-UHF result is in qualitatively aggreement with the FCI
one [see Fig. 4(b)], with the charge-carrier localization, however, being more
pronounced in the sS-UHF case. Again, it is remarkable that our FCI results
(including the crystal-field effect) for $\nu=3$ are in excellent agreement
with the experiment [36].

##  VI Summary

Taking into consideration (via a computational embedding scheme) both the
intra-moiré-QD charge-carrier interactions in the embedded MQD, as well as the
crystal-field effect from the surrounding moiré pockets (inter-moiré-QD
interactions), and using the full-configuration-interaction methodology, we
offered a complete interpretation of the counterintuitive experimental
observations reported in Ref. [36] in the context of moiré TMD superlattices
at integer fillings $\nu=2$ and 4. In particular, our results demonstrate that
these experimental observations reflect the interplay between two novel states
of matter; namely, (i) a supercrystal of sliding Wigner molecules 222When the
QD confinement possesses circular symmetry, the sliding WMs were earlier
referred to as rotating WMs or rotating electron molecules [14, 8]. for
unstrained moiré TMD materials (when the crystal field is commensurate with
the trilobal symmetry of the confining potential in each MQD) and (ii) a
supercrystal of pinned Wigner molecules when strain is involved and the
crystal field is non-commensurate with the trilobal symmetry of the confining
potential in each MQD. The case of $\nu=3$ is an exception, in that both
unstrained and strained cases produce a supercrystal of pinned WMs; this is
due to the congruence of intrinsic (that of the WM) and external (that of the
confining potential of the MQD) $C_{3}$ point-group symmetries.

Furthermore, we demonstrated that the UHF approach (invoked in Refs. [36, 56])
utterly fails to describe the genuinely quantum-mechanical supercrystal of
sliding WMs in the unstrained case, providing a qualitative agreement only in
the case of a supercrystal of pinned WMs.

The insights gained in this work through systematic comparisons, employing the
UHF method with exact diagonalization (FCI) calculations as benchmarks,
provide impetus and pathways for the development of effective strategies for
advancing economically feasible and accurate computational methodologies for
treating strongly-correlated many-body states in materials, including the
electron (hole) Wigner-molecular TMD superlattice systems studied here.

##  VII ackowledgements

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

##  Appendix A THE CONFIGURATION INTERACTION METHOD

The full configuration interaction (FCI) methodology has a long history,
starting in quantum chemistry; see Refs. [40, 45]. The method was adapted to
two dimensional problems and found extensive applications in the fields of
semiconductor quantum dots [8, 9, 42, 16, 14, 43, 26, 44] and of the
fractional quantum Hall effect [10, 22].

Our 2D FCI is described in our earlier publications. The reader will find a
comprehensive exposition in Appendix B of Ref. [26], where the method was
applied to GaAs double-quantum-dot quantum computer qubits. We specify that,
in the application to moiré DQDs, we keep similar space orbitals,
$\varphi_{j}(x,y)$, $j=1,2,\ldots,K$, that are employed in the building of the
single-particle basis of spin-orbitals used to construct the Slater
determinants $\Psi_{I}$, which span the many-body Hilbert space [see Eq. (B4)
in Ref. [26]; the index $I$ counts the Slater determinants]. The orbitals
$\varphi_{j}(x,y)$ are determined as solutions (in Cartesian coordinates) of
the auxiliary Hamiltonian

$\displaystyle H_{\rm aux}=\frac{{\bf
p}^{2}}{2m^{*}}+\frac{1}{2}m^{*}\omega_{0}^{2}(x^{2}+y^{2}),$ (6)

where $m^{*}\omega_{0}^{2}=16\pi^{2}v_{0}\cos(\phi)/a_{M}^{2}$, i.e., only the
isotropic parabolic (harmonic) contribution of the MQD confinement $V_{\rm
MQD}({\bf r})$ [see Eqs. (2) and (4) in the main text] is included.

The matrix elements
$\langle\varphi_{i}(x,y)|\sin(3\theta)r^{3}|\varphi_{j}(x,y)\rangle$ of the
anisotropic term in the MQD confinement are calculated analytically using the
algebraic language MATHEMATICA [57] and the Hermite-to-Laguerre (Cartesian-to-
polar) transformations listed in Ref. [58], whereas the matrix elements
$\langle\varphi_{i}(x,y)|H_{\rm CF}|\varphi_{j}(x,y)\rangle$ of the crystal
field are calculated numerically.

Following Ref. [26], we use a sparse-matrix eigensolver based on Implicitly
Restarted Arnoldi methods to diagonalize the many-body Hamiltonian in Eq. (3)
of the main text. We stress that the FCI conserves both the total spin and the
total-spin projection.

##  Appendix B THE SPIN-AND-SPACE UNRESTRICTED HARTREE-FOCK

Early on in the context of 2D materials, the spin-and-space unrestricted
Hartree-Fock (sS-UHF) was employed in Ref. [1] to describe formation of Wigner
molecules at the mean-field level. This methodology employs the Pople-Nesbet
equations [45, 14]. The sS-UHF WMs are self-consistent solutions of the Pople-
Nesbet equations that are obtained by relaxing both the total-spin and space
symmetry requirements. For a detailed description of the Pople-Nesbet
equations in the context of three-dimensional natural atoms and molecules, see
Ch. 3.8 in Ref. [45]. For a detailed description of the Pople-Nesbet equations
in the context of two-dimensional artificial atoms and semiconductor quantum
dots, see Sec. 2.1 of Ref. [14]. We note that the Pople-Nesbet equations
conserve the total-spin projection, but not the total spin. Convergence of the
self-consistent iterations was achieved in all cases by mixing the input and
output charge densities at each iteration step. The convergence criterion was
set to a difference of $10^{-12}$ meV between the input and output total UHF
energies at the same iteration step.

##  Appendix C CHARGE DENSITIES FROM FCI AND UHF WAVE FUNCTIONS

The FCI single-particle density (charge density) is the expectation value of a
one-body operator

$\rho({\bf r})=\langle\Phi^{\rm FCI}|\sum_{i=1}^{N}\delta({\bf r}-{\bf
r}_{i})|\Phi^{\rm FCI}\rangle,$ (7)

where $\Phi^{\rm FCI}$ denotes the many-body (multi-determinantal) FCI wave
function, namely,

$\Phi^{\rm FCI}({\bf r}_{1},\ldots,{\bf r}_{N})=\sum_{I}C_{I}\Psi_{I}({\bf
r}_{1},\ldots,{\bf r}_{N}),$ (8)

with $\Psi_{I}({\bf r)}$ denoting the Slater determinants that span the many-
body Hilbert space. The expansion coefficients $C_{I}$ are a byproduct of the
exact diagonalization of the many-body Hamiltonian $H_{\rm MB}$.

For the sS-UHF case, one substitutes $\Phi^{\rm FCI}$ in Eq. (7) with the
single-determinant, $\Psi^{\rm UHF}({\bf r})$, solution of the Pople-Nesbet
equations. $\Psi^{\rm UHF}({\bf r})$ is built out from the UHF spin-orbitals
whose space part has the form:

$\displaystyle u^{\alpha}_{i}=\sum_{\mu=1}^{K}{\cal C}^{\alpha}_{\mu
i}\varphi_{\mu},\;\;\;i=1,\ldots,K,$ (9)

and

$\displaystyle u^{\beta}_{i}=\sum_{\mu=1}^{K}{\cal C}^{\beta}_{\mu
i}\varphi_{\mu},\;\;\;i=1,\ldots,K,$ (10)

where the expansion coefficients ${\cal C}^{\alpha}_{\mu i}$ and ${\cal
C}^{\beta}_{\mu i}$ are solutions of the Pople-Nesbet equations, with $\alpha$
and $\beta$ denoting up and down spins, respectively.

The UHF expression for the charge density simplifies to:

$\displaystyle\rho^{\rm UHF}({\bf
r})=\sum_{i=1}^{N^{\alpha}}|u^{\alpha}_{i}({\bf
r})|^{2}+\sum_{j=1}^{N^{\beta}}|u^{\beta}_{j}({\bf r})|^{2},$ (11)

with ${N^{\alpha}}$, ${N^{\beta}}$ being the number of lowest-energy occupied
spin-up and spin-down orbitals, respectively.

## References

  * 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).
  * Egger  _et al._ [1999] R. Egger, W. Häusler, C. H. Mak, and H. Grabert, Crossover from Fermi liquid to Wigner molecule behavior in quantum dots, Phys. Rev. Lett. 82, 3320–3323 (1999).
  * Yannouleas and Landman [2000a] C. Yannouleas and U. Landman, Collective and Independent-Particle Motion in Two-Electron Artificial Atoms, Phys. Rev. Lett. 85, 1726–1729 (2000a).
  * Filinov  _et al._ [2001] A. V. Filinov, M. Bonitz, and Y. E. Lozovik, Wigner Crystallization in Mesoscopic 2D Electron Systems, Phys. Rev. Lett. 86, 3851–3854 (2001).
  * Yannouleas and Landman [2002a] C. Yannouleas and U. Landman, Strongly correlated wavefunctions for artificial atoms and molecules, Journal of Physics: Condensed Matter 14, L591–L598 (2002a).
  * Mikhailov [2002] S. A. Mikhailov, Two ground-state modifications of quantum-dot beryllium, Phys. Rev. B 66, 153313 (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).
  * 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 [2004] C. Yannouleas and U. Landman, Structural properties of electrons in quantum dots in high magnetic fields: Crystalline character of cusp states and excitation spectra, Phys. Rev. B 70, 235319 (2004).
  * Szafran  _et al._ [2004] B. Szafran, F. M. Peeters, S. Bednarek, and J. Adamowski, Anisotropic quantum dots: Correspondence between quantum and classical Wigner molecules, parity symmetry, and broken-symmetry states, Phys. Rev. B 69, 125344 (2004).
  * Romanovsky  _et al._ [2006] I. Romanovsky, C. Yannouleas, L. O. Baksmaty, and U. Landman, Bosonic molecules in rotating traps, Phys. Rev. Lett. 97, 090401 (2006).
  * Li  _et al._ [2006] Y. Li, C. Yannouleas, and U. Landman, From a few to many electrons in quantum dots under strong magnetic fields: Properties of rotating electron molecules with multiple rings, Phys. Rev. B 73, 075301 (2006).
  * 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).
  * Dai  _et al._ [2007] Z. Dai, J.-L. Zhu, N. Yang, and Y. Wang, Spin-dependent rotating Wigner molecules in quantum dots, Phys. Rev. B 76, 085308 (2007).
  * Li  _et al._ [2007] Y. Li, C. Yannouleas, and U. Landman, Three-electron anisotropic quantum dots in variable magnetic fields: Exact results for excitation spectra, spin structures, and entanglement, Phys. Rev. B 76, 245310 (2007).
  * Baksmaty  _et al._ [2007] L. O. Baksmaty, C. Yannouleas, and U. Landman, Rapidly rotating boson molecules with long- or short-range repulsion: An exact diagonalization study, Phys. Rev. A 75, 023620 (2007).
  * Ghosal  _et al._ [2007] A. Ghosal, A. D. Güçlü, C. J. Umrigar, D. Ullmo, and H. U. Baranger, Incipient Wigner localization in circular quantum dots, Phys. Rev. B 76, 085341 (2007).
  * Yang  _et al._ [2008] N. Yang, J.-L. Zhu, and Z. Dai, Rotating Wigner molecules and spin-related behaviors in quantum rings, Journal of Physics: Condensed Matter 20, 295202 (2008).
  * Romanovsky  _et al._ [2009] I. Romanovsky, C. Yannouleas, and U. Landman, Edge states in graphene quantum dots: Fractional quantum Hall effect analogies and differences at zero magnetic field, Phys. Rev. B 79, 075311 (2009).
  * Brandt  _et al._ [2015] B. B. Brandt, C. Yannouleas, and U. Landman, Double-well ultracold-fermions computational microscopy: Wave-function anatomy of attractive-pairing and Wigner-molecule entanglement and natural orbitals, Nano Letters 15, 7105–7111 (2015).
  * Yannouleas and Landman [2021] C. Yannouleas and U. Landman, Exact closed-form analytic wave functions in two dimensions: Contact-interacting fermionic spinful ultracold atoms in a rapidly rotating trap, Phys. Rev. Research 3, L032028 (2021).
  * Ercan  _et al._ [2021] H. E. Ercan, S. N. Coppersmith, and M. Friesen, Strong electron-electron interactions in Si/SiGe quantum dots, Phys. Rev. B 104, 235302 (2021).
  * Abadillo-Uriel  _et al._ [2021] J. C. Abadillo-Uriel, B. Martinez, M. Filippone, and Y.-M. Niquet, Two-body Wigner molecularization in asymmetric quantum dot spin qubits, Phys. Rev. B 104, 195305 (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).
  * 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).
  * Kalliakos  _et al._ [2008] S. Kalliakos, M. Rontani, V. Pellegrini, C. P. García, A. Pinczuk, G. Goldoni, E. Molinari, L. N. Pfeiffer, and K. W. West, A molecular state of correlated electrons in a quantum dot, Nature Physics 4, 467–471 (2008).
  * 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).
  * Jang  _et al._ [2023] W. Jang, J. Kim, J. Park, G. Kim, M.-K. Cho, H. Jang, S. Sim, B. Kang, H. Jung, V. Umansky, and D. Kim, Wigner-molecularization-enabled dynamic nuclear polarization, Nature Communications 14, 2948 (2023).
  * 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).
  * Palma  _et al._ [2023] F. D. Palma, F. Oppliger, W. Jang, S. Bosco, M. Janík, S. Calcaterra, G. Katsaros, G. Isella, D. Loss, and P. Scarlino, Strong hole-photon coupling in planar ge: probing the charge degree and wigner molecule states (2023), arXiv:2310.20661 [quant-ph] .
  * 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).
  * Li  _et al._ [2023] 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 (2023), arXiv:2312.07607 [cond-mat.mes-hall] .
  * Note [1] We note that the experimental results for the unstrained case in Ref. [36] were necessarily interpreted with formation of sliding WMs in isolated MQDs (for $N=2$ confined charge carriers, see Refs. [34, 36]; for $N=4$ confined charge carriers, see Ref. [34]), given the fact that the mean-field Hartree-Fock approach (that allows calculations for a large number of charge carriers distributed over several moiré pockets) is unreliable [35, 36] for describing sliding WMs.
  * Bethe [1929] H. Bethe, Termaufspaltung in Kristallen, Annalen der Physik 395, 133–208 (1929).
  * Van Vleck [1932] J. H. Van Vleck, Theory of the variations in paramagnetic anisotropy among different salts of the iron group, Phys. Rev. 41, 208–215 (1932).
  * Shavitt [1998] I. Shavitt, The history and evolution of configuration interaction, Molecular Physics 94, 3–17 (1998).
  * Yannouleas and Landman [2006] C. Yannouleas and U. Landman, Electron and boson clusters in confined geometries: Symmetry breaking in quantum dots and harmonic traps, Proceedings of the National Academy of Sciences 103, 10600–10605 (2006).
  * Rontani  _et al._ [2006] M. Rontani, C. Cavazzoni, D. Bellucci, and G. Goldoni, Full configuration interaction approach to the few-electron problem in artificial atoms, The Journal of Chemical Physics 124, 124102 (2006).
  * Li  _et al._ [2009] Y. Li, C. Yannouleas, and U. Landman, Artificial quantum-dot Helium molecules: Electronic spectra, spin structures, and Heisenberg clusters, Phys. Rev. B 80, 045326 (2009).
  * 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).
  * Szabo and Ostlund [1989] A. Szabo and N. S. Ostlund,  _Modern Quantum Chemistry_ (McGraw-Hill, New York, 1989) For the Slater-Condon rules, see Chap. 4.
  * Yannouleas and Landman [2002b] C. Yannouleas and U. Landman, Magnetic-field manipulation of chemical bonding in artificial molecules, International Journal of Quantum Chemistry 90, 699–708 (2002b).
  * Yannouleas and Landman [2000b] C. Yannouleas and U. Landman, Formation and control of electron molecules in artificial atoms: Impurity and magnetic-field effects, Phys. Rev. B 61, 15895–15904 (2000b).
  * Wu  _et al._ [2018] F. Wu, T. Lovorn, E. Tutuc, and A. H. MacDonald, Hubbard model physics in transition metal dichalcogenide moiré bands, Phys. Rev. Lett. 121, 026402 (2018).
  * Angeli and MacDonald [2021] M. Angeli and A. H. MacDonald, $\Gamma$-valley transition metal dichalcogenide moiré bands, Proceedings of the National Academy of Sciences 118, e2021826118 (2021).
  * Zhang  _et al._ [2020] Y. Zhang, N. F. Q. Yuan, and L. Fu, Moiré quantum chemistry: Charge transfer in transition metal dichalcogenide superlattices, Phys. Rev. B 102, 201115 (2020).
  * Wigner [1934] E. Wigner, On the interaction of electrons in metals, Phys. Rev. 46, 1002–1011 (1934).
  * Wigner [1938] E. Wigner, Effects of the electron interaction on the energy levels of electrons in metals, Trans. Faraday Soc. 34, 678–685 (1938).
  * Li  _et al._ [2021] H. Li, S. Li, E. C. Regan, D. Wang, W. Zhao, S. Kahn, K. Yumigeta, M. Blei, T. Taniguchi, K. Watanabe, S. Tongay, A. Zettl, M. F. Crommie, and F. Wang, Imaging two-dimensional generalized Wigner crystals, Nature 597, 650–654 (2021).
  * Tsui  _et al._ [2024] Y.-C. Tsui, M. He, Y. Hu, E. Lake, T. Wang, K. Watanabe, T. Taniguchi, M. P. Zaletel, and A. Yazdani, Direct observation of a magnetic-field-induced wigner crystal, Nature 628, 287–292 (2024).
  * Note [2] When the QD confinement possesses circular symmetry, the sliding WMs were earlier referred to as rotating WMs or rotating electron molecules [14, 8].
  * Reddy  _et al._ [2023] A. P. Reddy, T. Devakul, and L. Fu, Artificial Atoms, Wigner Molecules, and an Emergent Kagome Lattice in Semiconductor Moiré Superlattices, Phys. Rev. Lett. 131, 246501 (2023).
  * [57] Wolfram Research, Inc., Mathematica, Version 13.2, Champaign, IL, 2022\.
  * Kimel and Elias [1993] I. Kimel and L. Elias, Relations between Hermite and Laguerre Gaussian modes, IEEE Journal of Quantum Electronics 29, 2562–2567 (1993).