\section{Introduction}
Strong correlations are at the heart of many important phenomena in condensed matter systems, including unconventional superconductivity\cite{bednorz1986possible}, quantum magnetism\cite{fazekas1999lecture} and fractional quantum Hall states\cite{stormer1983quantized}. These phenomena have a wide range of applications in material design, energy science, and quantum information\cite{dagotto1994correlated}. The complexity of strongly correlated many-body systems does not allow us to apply traditional theoretical approaches based on perturbation theory, and requires using hard-core numerical techniques, including exact diagonalization\cite{lin1990exact}, quantum Monte Carlo\cite{foulkes2001quantum}, density-matrix renormalization group\cite{schollwock2005density}, \emph{etc}. However, these numerical techniques are limited to restricted conditions such as small size, high temperature, and low dimension. The pursuit for understanding strongly correlated systems in materials motivates new approaches that can overcome these restrictions.

In addition to conventional numerical techniques, analog quantum simulators offer a distinct solution. Specifically, cold-atom simulators in optical lattices have achieved great success in simulating interacting bosonic systems\cite{bloch2008many, lewenstein2012ultracold, gross2017quantum} and have recently begun exploring fermionic systems\cite{parsons2015site, greif2016site, parsons2016site, boll2016spin, brown2019bad,nichols2019spin}. Taking advantage of electrons as charged particles, solid-state quantum-dot simulators naturally incorporate the Coulomb interactions and provide an alternative for mimicking electronic many-body states in molecules and solids\cite{van2002electron, hanson2007spins, nielsen2007nanoscale, oguri2007kondo, thalineau2012few, seo2013charge, hensgens2017quantum, mukhopadhyay20182}. With the relatively easy accessibility of high orbitals and low temperatures, the quantum-dot simulators are promising to simulate a realistic system. Despite the experimental progress with these platforms for quantum simulation, the interpretation of the underlying physics is still at the stage of minimal models with estimated parameters\cite{hilker2017revealing, mazurenko2017cold, barthelemy2013quantum}. This limits the quantitative analysis of fine details of experiments and hinders extensions to more complicated models.

A solution to this problem might be readily available, if we turn to the fields of chemistry and material science, where atomic-basis-based \emph{ab initio} approaches have been well developed. The spirit of these approaches is the unbiased evaluation of all physical parameters from a given set of atomic ingredients. In the past half a century, \emph{ab initio} calculations have made great progress towards describing systems with increasing complexity. With the help of the Gaussian basis\cite{binkley1980self,ditchfield1971self,frisch1984self,hariharan1973influence,hehre1969self}, the computational cost has been largely reduced, making the simulation of large molecules possible. In addition to the basic Hartree-Fock method\cite{hartree1928wave,hartree1935self,fock1930naherungsmethode,fock1930selfconsistent,slater1930note}, many advanced post-Hartree-Fock wavefunction-based methods (coupled cluster\cite{shavitt2009many}, configuration interactions\cite{sherrill1999configuration,pople1987CI}, \emph{etc.}) and multi-reference methods\cite{hegarty1979application, eade1981direct, yamamoto1996direct, andersson1990second, andersson1992second} have been invented.  More recently, advanced computer architectures including graphical processing units (GPUs) have been widely exploited by quantum chemistry simulations, pushing the scale of calculation to even larger systems\cite{ufimtsev2008graphical,ufimtsev2008quantum, deprince2011coupled, hohenstein2015atomic, song2018reduced}. Though successful in chemistry, existing software packages are not compatible with quantum simulators: the state-of-the-art quantum chemistry calculations are based on existing atomic wavefunction bases; however, the tunability of quantum simulators requires the wavefunction basis being the eigenstates of given, arbitrary potential landscapes, which are obtained numerically during the calculation.

In small quantum-dot systems, initial progress has been made using fixed wavefunction bases. Early studies focused on the physics in a single parabolic quantum dot, whose basis wavefunctions are Fock-Darwin states -- a subset of Gaussian wavefunctions. With these analytical wavefunctions as bases, the many-body Hamiltonian can be easily computed and the ground-state solution can be obtained using exact diagonalization\cite{maksym1990quantum,merkt1991energy,pfannkuche1993comparison,hawrylak1993magnetoluminescence,yang1993addition,hawrylak1993single,palacios1994capacitance,wojs1995negatively,oaknin1995low,wojs1996charging,maksym1996eckardt,hawrylak1996magnetoexcitons,wojs1997theory,wojs1997spectral,eto1997electronic,maksym1998quantum,reimann2000formation,mikhailov2002quantum} and quantum Monte Carlo\cite{bolton1996fixed,harju1999many,harju1999wave,harju2002wigner,siljamaki2002various}. The simplified treatment was also extended to double-dot systems\cite{burkard1999coupled, li2009artificial,nielsen2012many,stepanenko2012singlet, yannouleas2016ultracold, brandt2017bottom}. Since the parabolic potentials cannot describe the ``crystal-field'' corrections -- the impact of neighboring dot potentials on the single-particle wavefunctions and site energies -- recent work has considered more realistic Gaussian potentials. In this case, density functional theory (DFT) and wavefunction based methods such as configuration interaction have been attempted, using a numerical wavefunction basis beyond the Fock-Darwin states\cite{abolfath2006real,stopa2006electronic,abolfath2006theory}. However, these are the largest quantum-dot systems that have been subject to \emph{ab initio} attempts. In trying to find a compromise between model accuracy and computational complexity, simulations of larger systems have been restricted to simpler toy models like the Hubbard and extended-Hubbard models\cite{stafford1994collective,stafford1997coherent,stafford1998coherent,kotlyar1998addition,kotlyar1998correlated}.

A recent experiment~\cite{dehollain2019nagaoka} showcased some of the power of quantum dot based simulators for studying quantum magnetism, by using a 2$\times$2 plaquette to investigate Nagaoka magnetism
-- magnetism induced by a single hole in a half-filled correlated electronic system. This phenomenon has been difficult to realize experimentally, in great part because of the correlated nature of the electronic system required to observe the physics of Nagaoka ferromagnetism\cite{nagaoka1966ferromagnetism}. The success of the experiment in Ref.~\onlinecite{dehollain2019nagaoka} relied on pushing the limits of the maximum achievable interaction strengths, as well as the minimum measurable energy gaps. The observed energy gap crucial for Nagaoka ferromagnetism is of the order of a few $\mu$eV in such a system, three orders of magnitude smaller than the level spacing between orbitals and the ground-state Coulomb interaction.

Given that Nagaoka ferromagnetism was proven in a single-band finite system, it is not obvious that this phenomenon should persist when the level spacing among different orbitals is well below the interaction scales, as is the case in the quantum-dot experiment by Dehollain \emph{et al}~\cite{dehollain2019nagaoka}. Thus, these system conditions require a precise numerical many-body approach in order to validate the experimental observations. Moreover, to reflect the tunability of quantum dots comparable with realistic experiments, the modeling with \emph{ab initio} inputs is also necessary.

For both of these purposes, we hereby introduce an \emph{ab initio} exact diagonalization framework to describe artificial quantum simulator systems consisting of multiple quantum dots. By calculating the wavefunctions in a given potential well and evaluating the one-center and two-center integrals, we construct the tight-binding Hamiltonian of the many-body system consisting of multiple interacting quantum dots. This calculation predicts the single-particle energies, along with various interaction energies, which are quantitatively consistent with experiments\cite{dehollain2019nagaoka}. Additionally, we applied the calculation on a plaquette system, reproducing the experimental conditions that led to the observation of the Nagaoka ferromagnetic ground state. The model again shows good agreement with the experimentally observed energy gaps, as well as with the observed robustness of the ferromagnetic state performed in the experiment\cite{dehollain2019nagaoka}.

The description of this model and calculation will gradually increase in complexity. In Sec.~\ref{sec:singleWell}, we first explain the single-well wavefunction basis and the numerical implementation that automatically generates the basis based on a given potential. After that, we present the derivation and implementation of many-body Hamiltonians in multiple quantum wells in Sec.~\ref{sec:manybodyHam}. By adjusting the model to represent a four-well system, in Sec.~\ref{sec:calculations} we then explore the quantum magnetism and especially the Nagaoka transition using the \emph{ab initio} exact diagonalization approach. Finally, we conclude and discuss the future directions of our approach in Sec.~\ref{sec:conclusion}.

\section{Single Electron in a Single-Quantum Well}\label{sec:singleWell}
\begin{figure}[!t]
\begin{center}
\includegraphics[width=\columnwidth]{Fig_singleWellSolution.jpg}
\caption{\label{fig:1} Solution of single-well wavefunctions for $V_0 = 100$ and $\Sigma =1$. (a) The Gaussian quantum well in two dimensions. (b) Eigen-energy solutions for all bound states in the quantum well of (a), with the colors denoting different angular quantum numbers. (c) Sample eigenstate wavefunctions for $(n,m)=(0,0)$, (4,0), (1,1), (4, 1), (1,5), and (2,7), respectively.
}
\end{center}
\end{figure}
To simulate the electrons trapped in a finite-width potential well, we consider a confining central potential with rotational symmetry. This confining potential mimics the combined impact of electrodes surrounding the quantum dot\cite{potentialshape}. Though a generic potential landscape, obtained by solving the Poisson equation, can be employed as an input in the calculation, we use the Gaussian potential $V(\rbf) = -V_0 e^{-{|\rbf|^2}/{2\sigma}}$ in this paper [see Fig.~\ref{fig:1}(a)], as a typical description of the finite-size quantum dot\cite{abolfath2006theory,abolfath2006real}. Here $\rbf = r\cos\phi\, \mathbf{e}_x + r\sin\phi\, \mathbf{e}_y$ is the spatial coordinate with respect to the center of quantum well. In contrast to an atomic potential, the quantum well has finite potential energy with no singularity; unlike the parabolic potential, the Gaussian potential has a finite width and finite number of bound states. The single-electron static \Schrdg\ equation is
\begin{equation}
\left[-\frac{\hbar^2}{2m_e^\star}\! \left(\!\frac{\partial^2}{\partial r^2} + \frac1{r}\frac{\partial}{\partial r} + \frac1{r^2}\frac{\partial^2}{\partial \phi^2}\!\right)\! +\! V(r)\right]\! \psi(\rbf) = E\psi(\rbf)\,,
\end{equation}
where $m_e^\star$ is the effective mass of electron in the two-dimensional electronic gas (2DEG). The equation can be simplified by separation of variables
\begin{equation}
    \psi(\rbf) = \frac{\chi(r)}{\sqrt{2\pi r}} e^{im\phi} = \frac{\chi(r)}{\sqrt{r}}\varphi(\phi),
\end{equation}
where the $\chi(r)$ and $\varphi(\phi)$ are the radius and angular wavefunctions. Denoting the radial quantum number as $n$ and angular quantum number as $m$, the set of $\{\chi_n(r)\}$ satisfies the normalization condition
\begin{equation}
    \int_0^\infty \chi_n(r)^* \chi_{n^\prime}(r) dr = \delta_{nn^\prime}\,.
\end{equation}
Then we obtain the radial differential equation
\begin{equation}
    -\frac{\hbar^2}{2m_e^\star}\left[\frac{d^2\chi}{dr^2} - \frac{m^2-1/4}{r^2} \chi\right] + V(r)\chi = E\chi\,.
\end{equation}
It can be numerically solved using the finite difference approximation.
Choosing the angular part being real for numerical convenience, we define the \textit{single-well} wavefunction as
\begin{equation}\label{eq:singleWellBasis}
    \psi_{nm}(\rbf) = \left\{\begin{array}{ll}
        \frac{\chi_n(r)}{\sqrt{2\pi r}}, & m=0 \\
        \frac{\chi_n(r)}{\sqrt{\pi r}} \cos(m\phi), & m>0\\
        \frac{\chi_n(r)}{\sqrt{\pi r}} \sin(m\phi), & m<0.
    \end{array}\right.
\end{equation}
These eigenstate wavefunctions define the 2D orbital $(n,m)$ quantum numbers, while the spin component will be introduced later. We label the single-well single-electron eigen-energy as $\varepsilon_{nm}$. As shown in Fig.~\ref{fig:1}(b), the energy levels are well separated near the ground state, but become denser at higher energies. This is typical in a finite potential well. Unlike a parabolic potential, there are finite number of bound states (denoted as $N_{\rm orbital}$) in a finite well.

The wavefunctions of the eigenstates also become more extended with the increase of energy, or equivalently quantum numbers. While $m$ determines the angular distribution of a wavefunction, $n$ gives the number of nodes along the radius. Figure~\ref{fig:1}(c) shows examples of a few eigenstate wavefunctions. The ground state $(n,m)=(0,0)$ is restricted to the center of the potential well with a Gaussian-like shape, while the high-energy states such as $(n,m)=(2,7)$ spread three times wider.

Different from 3D systems, the eigenstates of a 2D potential well have two-fold orbital degeneracy for all $|m|>0$ (\emph{i.e.},~$p$, $d$, $f$ orbitals in atomic notation). This degeneracy is maintained in a $C_4$ symmetric system. This rotational-symmetric shape of the potential well is a theoretical simplification. In reality, the confining potential is not perfectly symmetric and can deviate from the solution in Fig.~\ref{fig:1}, resulting in the level splitting of the degenerate states\cite{van2005spin}. However, as we will show in Sec.~\ref{sec:calculations}, the ideal model gives an adequate estimation of the experimentally measurable parameters, both qualitatively and quantitatively. This result indicates that single-well wavefunctions obtained from the rotational-symmetric potential well also form a good basis to expand local electronic states.

\section{Many-Body Model}\label{sec:manybodyHam}
With multiple potential wells, the general Hamiltonian for a many-body system among $N_{\rm well}$ wells is
\begin{equation}\label{eq:nonTighBindingHamiltonian}
\Ham = \sum_a^{N_e} \left[ -\frac{\hbar^2}{2m_e^\star} \nabla^2_a + \sum_i^{N_{\rm well}} V(\rbf_a - \textbf{R}_i)\right] + \sum_{a\neq b} \frac{e^2}{4\pi\epsilon|\rbf_a -\rbf_b|}\,,
\end{equation}
where the sum over $a$ and $b$ traverses the $N_e$ electrons, while the sum over $i$ traverses different potential wells. The first term is a sum with respect to each electron, which can be treated by separation of variables. Different from a chemistry problem, here the electrostatic potential $V(\rbf)$ is given by the electrodes and there is no need to introduce the Born-Oppenheimer approximation.

Following the linear combination of atomic orbitals (LCAO) approach in the electronic structure theory\cite{szabo2012modern}, we construct the basis using a superposition of the single-well wavefunctions
\begin{equation}\label{eq:orthTransform}
\tilde{\psi}_{\mu\sigma}(\rbf) = \sum_{\nu} X_{\nu\mu}\psi_{\nu\sigma}(\rbf)\,.
\end{equation}
For simplicity in notation, we collapse the coordinate and orbital indices as $\mu = (i,\alpha)$, and denote $\psi_{\mu\sigma}(\rbf) = \psi_\alpha(\rbf -\mathbf{R}_i)s(\sigma)$. The $s(\sigma)$ denotes the spin wavefunction which does not mix in the hybridization.
Since the single-well wavefunctions are truncated at a relatively high level, this linear combination does not span a complete spatial basis, but is enough for the ground-state calculation when the number of tracked orbitals is much larger than the number of occupied orbitals.
With the presence of multiple  wells, the single-well wavefunctions are no longer orthogonal. An orthonormalization should be applied in order to simplify the many-body \Schrdg\ equation. The overlap matrix among different single-well wavefunction basis is
\begin{equation}
\int d\rbf^3 \psi_{\mu\sigma}^*(\rbf )\psi_{\nu\sigma^\prime}(\rbf) = S_{\mu\nu}\delta_{\sigma\sigma^\prime}\,.
\end{equation}
Thus, the overlap matrix among the new basis functions is
\begin{eqnarray}
    \int d\rbf^3 \tilde\psi_{\mu_1\sigma_1}^*(\rbf)\tilde{\psi}_{\mu_2 \sigma_2}(\rbf)

    = X^\dagger SX\, \delta_{\sigma_1\sigma_2}\,.
\end{eqnarray}
By setting the requirement $X^\dagger SX= I$ and considering $S$ being positive-definitive, a standard choice\cite{szabo2012modern} is $X = S^{-1/2}$.
This selection results in a new orthonormal basis set $\{\tilde{\psi}_{\mu\sigma}(\rbf)\}$.

Representing the many-body wavefunction in the Fock space, spanned by the Slater determinants of $\{\tilde{\psi}_{\mu\sigma}(\rbf)\}$, we obtain the the second quantization of the many-body states\cite{helgaker2014molecular}:
\begin{eqnarray}
|\{i_k\alpha_k\sigma_k\}\rangle= c^\dagger_{i_N\alpha_N\sigma_N} \cdots c^\dagger_{i_2\alpha_2\sigma_2}c^\dagger_{i_1\alpha_1\sigma_1} |0\rangle\,,
\end{eqnarray}
where subscripted ``$i$''s denote the site indices labeling the quantum dots; ``$\alpha$''s denote the orbital indices $(n,m)$; ``$\sigma$''s denote the spin index.
To perform an exact diagonalization calculation, we further construct a second-quantized Hamiltonian, the generic form of which is
\begin{equation} \label{eq:tighBindingHamiltonian}
\Ham = \Ham_{\rm non-int} + \Ham_{\rm int}\,.
\end{equation}
The non-interacting part $\Ham_{\rm non-int}$ corresponds to the hopping of an electron across orbitals and wells, while the interacting part $\Ham_{\rm int}$ contains all the possible interactions between multiple electrons.
In the following subsections, we introduce the methodology and approximations to evaluate these two parts of Hamiltonians.

\subsection{Non-Interacting Part of Hamiltonian}
The non-interacting part of the Hamiltonian is quadratic in fermionic operators
\begin{equation} \label{eq:noninteractingHamiltonian}
\Ham_{\rm non-int} = \sum_{i\alpha\sigma}\sum_{j\beta\sigma} t_{ij\atop \alpha\beta} c_{i\alpha\sigma}^\dagger c_{j\beta\sigma}\,.
\end{equation}
For any Fock state, this quadratic term can be evaluated using single-electron states. To simplify the calculation, we can first evaluate the matrix elements using the original non-orthogonal basis $\{\psi_{\mu\sigma}(r)\}$, obtaining a matrix $h_{\mu\nu}$, which we can transform into the orthonormal basis. This results in
\begin{eqnarray}\label{eq:atomicOrbitalTTerms}
h_{\mu\nu}&=& \int d\rbf^3 \psi_{\mu\sigma}^*(\rbf)\left[ -\frac{\hbar^2}{2m_e^\star} \nabla^2 + \sum_i V(\rbf - \textbf{R}_i)\right]\psi_{\nu\sigma}(\rbf)\nonumber\\
&=& \varepsilon_\nu + \int d\rbf^3 \psi_{\mu\sigma}^*(\rbf)\sum_{j\neq i_\nu} V(\rbf - \textbf{R}_j)\psi_{\nu\sigma}(\rbf)\,,
\end{eqnarray}
for any spin $\sigma$.
The diagonal terms of $h_{\mu\nu}$ define the site energies associated with each (single-well) orbital $\nu$. Note, that this energy is not equal to the bare eigen-energy $\varepsilon_\nu$ in a single well, since the second term also has a finite diagonal contribution. This is an analog of the ``crystal field''. The off-diagonal terms in $h_{\mu\nu}$ define the hybridization between different orbitals.

The transformation into orthonormal basis is done by substituting Eq.~\eqref{eq:orthTransform} into \eqref{eq:atomicOrbitalTTerms}, resulting in
\begin{eqnarray}\label{eq:defTTerms}
t_{\mu\nu}&=&\langle i\alpha\sigma |\Ham_{\rm non-int} |j\beta\sigma\rangle\nonumber\\
 &=& \sum_{\mu^\prime\nu^\prime}\int d\rbf^3 X^*_{\mu^\prime\mu}\psi_{\mu^\prime\sigma}(\rbf)\Ham
 \psi_{\nu^\prime\sigma}(\rbf)X_{\nu^\prime\nu}\,
\end{eqnarray}
Here $t_{\mu\nu}$ defines the site energy (diagonal) and hybridization (off-diagonal) of the orthonormal orbitals, which appears in Eq.~\eqref{eq:tighBindingHamiltonian}. Due to the impact of the ``crystal field'' and hybridization, the energy distribution of a multi-well system can be dramatically different from the single-well solution.

\begin{figure}[!b]
\begin{center}
\includegraphics[width=\columnwidth]{Fig_excludedTermsCartoon.jpg}
\caption{\label{fig:multipletexcludedTerms} The interaction terms that are ignored in the tight-binding Hamiltonian: (a) the density-dependent hopping and (b) the scattering terms involving more than two orbitals with different site energies.}
\end{center}
\end{figure}

\subsection{Interacting Part of Hamiltonian}\label{sec:interactiongHam}
As a typical choice in condensed matter, we restrict the interaction part $\Ham_{\rm int}$ to four-fermion terms\cite{szabo2012modern}. As specified in Eq.~\eqref{eq:genericFourFermionTerms}, the generic second-quantized four-fermion term contains an enumeration of four coordinate indices ($i_1,i_2,j_1,j_2$), four orbital indices ($\alpha_1, \alpha_2, \beta_1, \beta_2$) and four spin indices ($\sigma_1, \sigma_2, \sigma^\prime_1,\sigma^\prime_2$). This results in $16N_{\rm well}^4 N_{\rm orbital}^4$ interaction terms, whose general expression is shown in Eq.~\eqref{eq:generalVterms}. The bottleneck of the computation is the evaluation of the interaction parameters by numerical integration. Therefore, the setup of a model involving all combinatorial possibilities is currently beyond our capability. Hence, we introduce several common approximations to reduce the number of independent variables. Firstly, without relativistic effects, the Coulomb interaction is independent of spin; therefore, $\{\sigma^\prime_1,\sigma^\prime_2\} = \{\sigma_1,\sigma_2\}$. Secondly, due to the two-body nature of the interaction, one-center and two-center integrals dominate, whereas terms with more centers decay exponentially for well-separated wells. Dropping these multi-center terms implies the assumption that the geometric coordinates \{$i_1$, $i_2$, $j_1$, $j_2$\} can take at most two values. Thirdly, to further reduce the complexity, we restrict the interaction terms to ``perfectly'' resonant processes\cite{resonance}, which strictly speaking is only fully justified when level splitting is much larger than the interaction energy scales. For example, we neglect two generic classes of interactions: the density-dependent hopping and the scattering terms involving more than two orbitals [see Fig.~\ref{fig:multipletexcludedTerms}]. These terms are important in some cold-atom systems where individual energy scales are controllable, but become non-resonant in our model due to the strong interaction and unequal spacing between energy levels\cite{nonresonance}. The omission of these two non-resonant processes is equivalent to setting $\{i_1,i_2\}=\{j_1,j_2\}$ and restricting each four-fermion interaction term to at most two orbital indices. The above approximations significantly reduce the complexity of the model and have been a common strategy in solid state\cite{dagotto2013nanoscale}.

\begin{figure}[!t]
\begin{center}
\includegraphics[width=\columnwidth]{Fig_onSiteMultipletCartoon.jpg}
\caption{\label{fig:multipletOnsite} On-site interactions within one quantum dot: Hubbard $U$, inter-orbital Hubbard $U^\prime$ (and its spin-anti-parallel form $\bar{U}^\prime$), and Hund's exchange $J$ (and its spin-anti-parallel form $\bar{J}$).}
\end{center}
\end{figure}

After these simplifications, the interacting part of the Hamiltonian can be decomposed as the on-site and (two-center) long-range parts
\begin{equation} \label{eq:decomposingIntHamiltonian}
\Ham_{\rm int} = \sum_i \Ham_i^{\rm(OS)} + \sum_{ij} \Ham_{ij}^{\rm(LR)}.
\end{equation}
The standard derivations of each term in the interacting Hamiltonian is present in Appendix~\ref{app:interactiongHam}.
The on-site interaction Hamiltonian can be written as
\begin{eqnarray}\label{eq:multipletHam}
\Ham^{\rm (OS)}_i\!&=&\!\frac12\sum_{\alpha\sigma} U_{\alpha} n_{\alpha\bar{\sigma}}n_{\alpha\sigma} + \frac12\!\sum_{\alpha_1\neq\alpha_2}\!\sum_{\sigma_1,\sigma_2} U^\prime_{\alpha_1 \alpha_2} n_{\alpha_2\sigma_2}n_{\alpha_1\sigma_1} \nonumber\\
&& +  \frac12\sum_{\alpha_1\neq\alpha_2}\sum_{\sigma_1,\sigma_2} J_{\alpha_1 \alpha_2}
 c^\dagger_{\alpha_2\sigma_1}c^\dagger_{\alpha_1\sigma_2}  c_{\alpha_2\sigma_2}c_{\alpha_1\sigma_1} \,.
\end{eqnarray}
This is the known as the \textit{multiplet model}, widely used to describe the valence electrons in the transition metal systems\cite{dagotto2013nanoscale}.
The corresponding scattering processes are sketched in Fig.~\ref{fig:multipletOnsite}. For convenience, the site index is removed on the right-hand side, while in an inhomogeneous system [such as the modulations in Sec.~\ref{sec:calculations}], one should consider it specifically for each individual site.
Due to a symmetry consideration [see discussions in Appendix~\ref{app:integral}], it is usually convenient to calculate the interaction parameters using the original single-well basis obtained from Eq.~\eqref{eq:singleWellBasis}, through
\begin{equation}\label{eq:expXi}
    \Xi_{\mu\!_1\nu\!_1\atop\mu\!_2\nu\!_2}\!=\! \iint\! d\rbf_1^dd\rbf_2^d\, W(|\rbf_1\! -\!\rbf_2|)\psi_{\mu_1}\!(\!\rbf_1\!)^*\psi_{\mu_2}\!(\!\rbf_2\!)^*\psi_{\nu_1}\!(\!\rbf_1\!)\psi_{\nu_2}\!(\!\rbf_2\!)\,.
\end{equation}
Here $W(|\rbf_1\! -\!\rbf_2|)={e^2}/{4\pi\epsilon|\rbf_1 -\rbf_2|}$ is the two-electron Coulomb repulsion.
Note, that here we have taken the compact notation $\mu = (j,\beta)$ introduced above and have omitted the spin indices as they do not affect the spatial integral. Then using Eq.~\eqref{eq:orthTransform}, we have
\begin{equation}\label{eq:onsiteParams}\begin{aligned}
    U^\prime_{\alpha_1 \alpha_2} &=\sum_{\mu_1,\mu_2}\sum_{\nu_1,\nu_2}X_{\mu_1a_1}^*X_{\nu_1 a_1}X_{\mu_2a_2}^*X_{\nu_2a_2}\Xi_{\mu_1\nu_1\atop\mu_2\nu_2}\\
    J_{\alpha_1 \alpha_2} &= \sum_{\mu_1,\mu_2}\sum_{\nu_1,\nu_2}X_{\mu_1a_2}^*X_{\nu_1 a_1}X_{\mu_2a_1}^*X_{\nu_2a_2}\Xi_{\mu_1\nu_1\atop\mu_2\nu_2}.
\end{aligned} \end{equation}
These parameters define the on-site interactions among the orthonormal orbitals.

\begin{figure}[!t]
\begin{center}
\includegraphics[width=0.95\columnwidth]{Fig_offSiteMultipletCartoon.jpg}
\caption{\label{fig:multipletOffsite} Long-range interactions between two quantum dots: direct Coulomb interaction $V$, long-range Hund's exchange $K$, correlated on-site exchange $V^\prime$ and correlated off-site exchange $K^\prime$.
}
\end{center}
\end{figure}

Similarly, the long-range interactions are written as
\begin{eqnarray}\label{eq:longrangeHam}
\Ham^{\rm (LR)}_{ij}&=&\frac12\sum_{\alpha\sigma}\sum_{\beta\sigma^\prime} V_{\alpha\beta} n_{i\alpha\sigma}n_{j\beta\sigma^\prime}\nonumber\\
&& + \frac12\sum_{\alpha\beta}\sum_{\sigma\sigma^\prime} K_{\alpha\beta}
 c^\dagger_{j\beta\sigma}c^\dagger_{i\alpha\sigma^\prime}  c_{j\beta\sigma^\prime}c_{i\alpha\sigma}\nonumber\\
 &&+ \frac12\sum_{\alpha\neq\beta}\sum_{\sigma\sigma^\prime} V^\prime_{\alpha\beta}
 c^\dagger_{i\beta\sigma}c^\dagger_{j\alpha\sigma^\prime}  c_{j\beta\sigma^\prime}c_{i\alpha\sigma}\nonumber\\
 && + \frac12\sum_{\alpha\neq\beta}\sum_{\sigma\sigma^\prime} K^\prime_{\alpha\beta}
 c^\dagger_{j\alpha\sigma}c^\dagger_{i\beta\sigma^\prime}  c_{j\beta\sigma^\prime}c_{i\alpha\sigma}\nonumber\\
 && + \frac12\sum_{\alpha\neq\beta}\sum_{\sigma\neq\sigma^\prime} K^{\prime\prime}_{\alpha\beta}
 c^\dagger_{i\beta\sigma}c^\dagger_{j\beta\sigma^\prime}  c_{j\alpha\sigma^\prime}c_{i\alpha\sigma}\,.
\end{eqnarray}
As sketched in Fig.~\ref{fig:multipletOffsite}, $V_{\alpha\beta}$ represents a direct Coulomb interaction and $K_{\alpha\beta}$ is the corresponding exchange interaction; similarly, $V^\prime_{\alpha\beta}$ is the correlation between two on-site exchange interactions, while $K^\prime_{\alpha\beta}$ is the correlation between off-site exchange. The $K^{\prime\prime}_{\alpha\beta}$ term is an analog of the pair-hopping term and is also ignored here.
The expressions for the relevant long-range terms are
\begin{eqnarray}\label{eq:longrangeParams}
    V_{\alpha\beta}^{ij} &=& \sum_{\mu_1,\mu_2}\sum_{\nu_1,\nu_2}X_{\mu_1(i\alpha)}^*X_{\nu_1(i\alpha)}X_{\mu_2(j\beta)}^*X_{\nu_2(j\beta)}\Xi_{\mu_1\nu_1\atop\mu_2\nu_2}\nonumber\\
    K_{\alpha\beta}^{ij} &=& \sum_{\mu_1,\mu_2}\sum_{\nu_1,\nu_2}X_{\mu_1(j\beta)}^*X_{\nu_1 (i\alpha)}X_{\mu_2(i\alpha)}^*X_{\nu_2(j\beta)}\Xi_{\mu_1\nu_1\atop\mu_2\nu_2}\nonumber\\
    V_{\alpha\beta}^{ij\prime} &=& \sum_{\mu_1,\mu_2}\sum_{\nu_1,\nu_2}X_{\mu_1(i\beta)}^*X_{\nu_1(i\alpha)}X_{\mu_2(j\alpha)}^*X_{\nu_2(j\beta)}\Xi_{\mu_1\nu_1\atop\mu_2\nu_2}\nonumber\\
    K_{\alpha\beta}^{ij\prime} &=& \sum_{\mu_1,\mu_2}\sum_{\nu_1,\nu_2}X_{\mu_1(j\alpha)}^*X_{\nu_1 (i\alpha)}X_{\mu_2(i\beta)}^*X_{\nu_2(j\beta)}\Xi_{\mu_1\nu_1\atop\mu_2\nu_2}\,.
\end{eqnarray}
Note, the long-range interaction has contributions from both direct long-range integrals (for two-center $\mu_i$ and $\nu_i$ indices), and indirect hybridized on-site integrals (for one-center $\mu_i$ and $\nu_i$ indices). With well-separated quantum dots, the long-range interactions are typically much smaller than the on-site interactions. That being said, $V\ll U$, $K\ll J$, and $V^\prime$ and $K^\prime$ are even smaller compared to $V$ and $K$. Due to the orbital match of on-site wavefunctions, the $V$ terms are expected to dominate in the long-range interactions. However, for the study of Nagaoka ferromagnetism in the plaquette (see Sec. ~\ref{sec:calculations}), it is necessary to consider all of these long-range parameters, since the effects we want to observe can depend significantly on the superfine structures.

\begin{figure}[!t]
\begin{center}
\includegraphics[width=\columnwidth]{Fig_multipletCartoon.jpg}
\caption{\label{fig:multiplet} Cartoon for multi-orbital Nagaoka transition in a four-dot system. For moderate effective interaction, the multiplets in each quantum well form an overall low-spin state, with total spin $S=1/2$ (left). In contrast, a large interaction relative to the tunneling gives a Nagaoka FM state (right). The shaded surfaces denote the potential wells, while the white dots denote the single-well energy levels, which are slightly shifted according to different angular quantum numbers. The spin configuration is a conceptual sketch instead of a realistic solution.
}
\end{center}
\end{figure}
\section{Simulation of Four-Well Quantum Dot System: Probing Nagaoka Magnetism}\label{sec:calculations}
The explicit expressions for the tight-binding parameters described above allow one to fully diagonalize many-body electronic systems with multiple quantum dots. We will use this methodology to investigate the physics described by Nagaoka\cite{nagaoka1966ferromagnetism}, applied to a multi-orbital, 2$\times$2 system. Specifically, we study a system with three electrons in a four-site plaquette, which realizes the condition of a single hole in a Mott insulator where for a single orbital per site Nagaoka proved that the ground state must be ferromagnetic in the limit of large interaction strength. As sketched in Fig.~\ref{fig:multiplet}, with the total electron occupation less than the number of quantum dots, the multiplets on each quantum dot interact with each other and are expected to yield an effective collective spin configuration. If the multi-orbital system has similar behaviors to those of a single-band system formed by those multiplets, we expect it to display a high-spin-low-spin transition at various model parameters: with large enough interaction relative to the tunneling, we expect the Nagaoka mechanism to yield a ferromagnetic (FM) high-spin ground state; however, with moderate interactions, the system becomes a doped Mott insulator with a low-spin ground-state configuration, which corresponds to an anti-ferromagnetic state in the thermodynamic limit\cite{von2010probing}.

A recent experiment has studied Nagaoka magnetism using a quantum dot array in a $2\times2$ plaquette configuration\cite{dehollain2019nagaoka}. For a great part of the analysis in that work, a single-band extended Hubbard model with fitted parameters was used to model the system, obtaining results that seem to describe most of the experimental observations accurately. However, the fact that the experimentally observed level spacing between the two lowest orbitals is smaller than the electronic interaction raises the question of whether the system is adequately described by the single-band model. In this section, we use the \emph{ab initio} exact diagonalization approach described above to extract the precise many-body model of the $2\times2$ quantum dot plaquette and quantitatively reproduce the Nagaoka conditions that were explored with the experimental system.

\subsection{Evaluation of Model Parameters}
To compare with a realistic system, we first discuss the typical values of parameters. The gate-electrode structure of the experimental device was lithographically designed to define quantum dot wells on the scale of 100\,nm\cite{mukhopadhyay20182,dehollain2019nagaoka}. Therefore, we set our spatial units of the lattice constant $a_0=100$\,nm and Gaussian potential width $\Sigma = 100$\,nm. Considering the effective mass of electrons in a GaAs/AlGaAs  2DEG is $m_e^\star \approx 0.067\,m_e$, the natural energy unit corresponds to ${\hbar^2}/{a_0^2m_e^\star}\approx 0.114$\,meV. Applying this scale to the eigen-spectrum solved in Fig.~\ref{fig:1} ({i.e.},~$V_0 = 100\, {\hbar^2}/{a_0^2m_e^\star} = 11.4$\,meV), we obtain the ground state to first excited state level spacing $\Delta E = \varepsilon_1-\varepsilon_0 \approx 0.75$\,meV, comparable to the experimental observation of $\sim 1$\,meV.

\begin{figure}[!b]
\begin{center}
\includegraphics[width=\columnwidth]{Fig_levelSpaceandU_vs_shape.jpg}
\caption{\label{fig:interaction} (a) Ground-first-excited state level spacing, (b) ground-state interaction $U$ and (c) effective interaction as a function of the depth $V_0$ and the width $\Sigma$ of quantum well. The calculation is obtained on a single quantum well without hybridization.
}
\end{center}
\end{figure}

The evaluation of the electron-electron interaction requires a specific value of the dielectric constant $\epsilon$, which is ideally 12.9 in GaAs. It is known that the presence of metallic gate electrodes in the vicinity of the 2DEG has the effect of increasing $\epsilon$. However, the precise evaluation of $\epsilon$ is challenging. Instead, we rely on the value of the addition energy, which has been accurately estimated by experiments to be 2.9 meV, and select an $\epsilon$ that results in reasonable interaction values. Taking $\epsilon = 20$ into the solution of $V_0 = 11.4$\,meV mentioned above gives the ground-state $U_0\approx2.34$\,meV and the ground-excited-state $U^\prime_{01} \approx 1.92$\,meV. Note, these are the intrinsic model parameters in the many-body Hamiltonian. A typical experimental estimation of this Hubbard interaction is obtained by measuring the addition energy. Due to the orbital mixture when $\Delta E < U$ and the fact that excited-state wavefunctions are spatially wider, the experimentally measured ``effective interaction'' strength is slightly smaller than the model parameters $U$ and $U^\prime$. Figure~\ref{fig:interaction} gives an example of level spacing $\Delta E$, ground-state Hubbard $U$ and effective interaction calculated in a single-well system with different shape parameters.

The long-range interactions are much smaller than the on-site ones. Specifically for $d=210$\,nm, the long-range interaction $V$ ranges from 0.22\,meV to 0.4\,meV depending on the orbitals; $K$ and $V^\prime$ are on the order of or below 1\,$\mu$eV; the off-site exchange correlation $K^\prime$s are even lower, on the order of 0.1 or 0.01\,$\mu$eV. These terms form higher-order corrections to the multiplet model of Eq.~\eqref{eq:multipletHam}. As shown in Table \ref{tab:modelComp}, only the long-range Coulomb interaction $V$ obviously affects the ground-state energy, by order of 1 meV, while others contribute to $\sim0.01$\,meV. However, as stated before and now made clear in Table \ref{tab:modelComp}, the strong interaction condition results in a high-spin to low-spin state energy gap--which we refer to as the Nagaoka gap, that is on the scale of $\mu$eV. The precise value of the Nagaoka gap depends on the details of the microscopic parameters such as the confining potential for electrons and the many-body interactions. Therefore, every long-range term provides a non-negligible contribution to the Nagaoka gap. Noticeably, the $K$ terms have larger contributions to the Nagaoka gap than $V$, although it is inconsistent to include only one of them because it is the combination of both that obeys the exchange relation in Eq.~\eqref{eq:generalVterms}. A closer inspection of the dependence of the Nagaoka gap size on various models -- in particular the contrast between $t$-$U$-$J$ and $t$-$U$-$J$-$V$-$K$-$V^\prime$-$K^\prime$ -- indicates that the long-range Hund's exchange only contributes $\sim 23\%$ of the ferromagnetic effect, with the Nagaoka mechanism dominating. Distinguishing these two contributions is only possible in a multi-band model. This quantitative assignment gives further confirmation that the experimental result in Ref.~\onlinecite{dehollain2019nagaoka} is indeed caused by a Nagaoka-like mechanism.

\begin{table}
    \begin{tabular}{c|c|c}
    \hline \hline
 Model& Ground-State Energy  & Nagaoka Gap  \\
 \hline
 $t$-$U$-$J$  & -43.579950\,meV & 2.213\,$\mu$eV \\
 $t$-$U$-$J$-$V$ & -42.576572\,meV & 2.318\,$\mu$eV\\
 $t$-$U$-$J$-$V$-$K$ & -42.558866\,meV & 2.775\,$\mu$eV \\
 $t$-$U$-$J$-$V$-$K$-$V^\prime$-$K^\prime$ & -42.558912\,meV & 2.868\,$\mu$eV\\
 \hline  \hline
\end{tabular}
\caption{\label{tab:modelComp} Effect of system parameters (definition of these parameters can be found in the Sec.~\ref{sec:interactiongHam}) on ground-state energies and the Nagaoka gap obtained by various models for $d=210$\,nm. The calculations are performed on a four-dot system with three electrons, and the ground states of all models listed here are high-spin states.}
\end{table}

\begin{figure}[!t]
\begin{center}
\includegraphics[width=\columnwidth]{Fig_effectiveHopping.jpg}
\caption{\label{fig:bandwidth} The effective hopping $t$ estimated by a quarter of the single-particle bandwidth calculated for various distances in a $2\times 2$ plaquette.
}
\end{center}
\end{figure}
The hybridizations, or tunneling terms, vary among different orbitals and are exponentially dependent on the distance between quantum dot potential wells. Since the single-well ground-state wavefunctions are most localized, the hybridizations between neighboring-quantum-dot (single-well) ground states are extremely small ($\sim0.06\,\mu$eV for $d=210$\,nm). However, these local orbitals and tunnelings among them $h_{\mu\nu}$ are non-physical: they are nothing but mathematical tools to solve the many-body problem\cite{coherence}. In reality, the ``crystal field'' and wavefunction orthogonalization cause heavy hybridization between the (single-well) ground state and excited states -- the maximum of which can be close to $\Delta E$. These high-level excited states can contribute a $\sim 0.5$\,meV hopping amplitude between neighboring quantum wells. Therefore, the experimentally measurable \emph{effective tunneling} across low-energy states is the result of a superposition of all different conceptual paths.

Following this philosophy, the effective hopping $t$ can be simply extracted from the single-particle bandwidth in the entire multi-well system. If we only consider nearest-neighbor tunneling, the low-energy band structure of a $2\times 2$ plaquette takes the form $E(\theta) = -2t\cos \theta$ for $\theta = \{0, \pi/2, \pi, 3\pi/2\}$. Therefore, the width of the lowest band (the lowest four states) in a single-electron system gives an estimation of $4t$. Figure~\ref{fig:bandwidth} shows the extracted values of $t$ for different neighboring-dot distances. In the experimental device, the inter-dot tunneling can be tuned to the range of 0-40\,$\mu$eV\cite{mukhopadhyay20182, dehollain2019nagaoka}, which in the \emph{ab initio} model corresponds to a range of distances $d=210-240$\,nm. This is fairly consistent with the lithographically designed inter-dot distance of 150\,nm, which is also an approximation, since the actual inter-dot distance in the experiments is not measurable.

We emphasize that the above model parameters (including the first excited-state level spacing $\Delta E$, the ground-state and ground-excited-state Coulomb interaction $U_0$ and $U_{01}^\prime$, long-range Coulomb interaction $V$, the effective tunneling $t$, and the Nagaoka gap $\Delta$) evaluated from our \emph{ab initio} calculation using only very limited experimental input  match quantitatively with the experiment in Ref.~\onlinecite{dehollain2019nagaoka}. Therefore, we believe the \emph{ab initio} calculation serves the purpose of predicting model parameters in a quantitative level based only on given potential landscapes.

To simulate the correlated Nagaoka physics in multiple quantum dots, we perform the calculation in a microcanonical ensemble, with three electrons in a four-well system, and focus on the ground-state properties. The evaluation of single-well eigenstates and the integration are performed on a grid with a spacing of 1\,nm. To simplify the calculation, we keep 15 orbitals in each quantum dot, which span a $\sim 5$\,meV energy range. As this range is much larger than both $U$ and $t$, we believe that the level mixture above this truncation can be ignored\cite{convergence}. We perform exact diagonalization to solve this 60-orbital spinful system, using the parallel Arnoldi approach\cite{lehoucq1998arpack, jia2017paradeisos}.

\begin{figure}[!t]
\begin{center}
\includegraphics[width=\columnwidth]{Fig_distanceDep.jpg}
\caption{\label{fig:distDep} (a) The Nagaoka gap and (b) the ground-state energy of three electrons in four quantum dots, as a function of the distance $d$. The red open circles denote the low-spin ground states, while the blue dots denote the high-spin ground states. The size of the data points reflects the energy difference between the lowest low-spin and high-spin states in a logarithmic scale.
}
\end{center}
\end{figure}
\subsection{Distance Dependence}
Having selected the quantum dot potential well parameters, we first study the ground state properties as a function of the distance between neighboring dots in the plaquette. As shown in Fig.~\ref{fig:distDep}(b), the energy increases monotonically when the quantum dot separation is increased from 200\,nm to 280\,nm. This is a consequence of the crystal field renormalization of the site energies. As the dot separation becomes large enough to make the long-range interactions negligible, the electrons can no longer lower their energy by delocalizing, and the ground-state energy saturates towards $\sim 30$\,meV. This energy corresponds to each of the three electrons occupying the ground state of a quantum well independently.

Interestingly, the ground-state configuration switches from a high- to a low-spin state at $d\gtrsim 206$\,nm. This is a feature of the Nagaoka effect applied to finite-size lattices, which have access to regimes outside of the thermodynamic limit ($U/t \rightarrow \infty$) where Nagaoka made the original prediction. Increasing the distance between dots effectively suppresses $t$ and long-range interactions, but changes little of the on-site interactions. At small enough effective tunneling with large enough distance, the $U\gg t$ condition is reached at some point. Such a Nagaoka effect was originally predicted for a single hole in a half-filled Hubbard model in the thermodynamic limit, where the transition occurs at an infinite $U/t$ ratio. However, this critical ratio becomes finite for a finite cluster, since the underlying physics reflected by the Nagaoka transition is a $t$ versus $N_{\rm well}\times J$ competition. This phenomenon was previously shown (and proven) in a single-band Hubbard\cite{nagaoka1966ferromagnetism, mattis2003eigenvalues} and extended Hubbard models\cite{kollar1996ferromagnetism}. Here we show its validity in a multi-orbital system.

\begin{figure}[!b]
\begin{center}
\includegraphics[width=\columnwidth]{Fig_detuneGroundState_new.jpg}
\caption{\label{fig:detuningGroundState} (a) Schematic illustrating the potential detuning applied on one of the four quantum dot potential wells. (b) The ground-state energy for the entire system as a function of the potential detuning $dV$, calculated for various distances $d$. The gray lines denote the energy drop with slope 1.
}
\end{center}
\end{figure}

As shown in Fig.~\ref{fig:distDep}(a), the Nagaoka gap switches to positive at $d> 206$\,nm and reaches a maximum at $d\sim 210$\,nm. With larger distances, the Nagaoka gap starts to decrease as the correlations among electrons in different wells diminish. We select $d= 210$\,nm as the default geometric setup for the following calculations. In this case, the absolute value of the Nagaoka gap is 2.87\,$\mu$eV, consistent with an estimation in Ref.~\onlinecite{dehollain2019nagaoka} through a comparison between experimentally measured parameters with a fitted single-band model.

\subsection{Potential Detuning}

In addition to investigating the Nagaoka transition as a function of separation between the dots, we demonstrate that the low-spin-high-spin transition can also be driven by varying the potential of a single well, which reflects the robustness of the magnetism against disorder. As shown in Fig.~\ref{fig:detuningGroundState}(a), we vary the depth $V_0$ of one of the wells by a positive or negative $dV$, which results in unbalanced site-energies. More broadly, the change of all eigenstates associated with this particular well affects the hybridization and interaction parameters. These changes are all captured in the \emph{ab initio} calculation.

The results from this study, shown in Fig.~\ref{fig:detuningGroundState}(b), give some expected, but also some unexpected outcomes. A first observation is that the total energy of the system is lowered as the selected well is made deeper, and the Nagaoka condition breaks when the well becomes sufficiently shallow or deep. Surprisingly though, the slope of such energy decrease varies when $dV$ switches from positive to negative. Additionally, there is an asymmetry in the robustness of the Nagaoka state, between positive and negative detuning, which was also observed in the experiment\cite{dehollain2019nagaoka}. Taking the $d=210$\,nm system as an example, at $dV=0$, the ground state is the Nagaoka high-spin state discussed above; when the potential detuning is $dV\!=\!0.11$\,meV or $dV\!=\!-0.07$\,meV, the system undergoes a transition to the low-spin ground state. The asymmetric behaviors indicate that the transitions at positive and negative $dV$s have a different nature.

For $dV>0$, the detuned well is deeper, lowering the energy barrier for a doubly-occupied state (sometimes called doublon) and accordingly increasing the spin-exchange energy $J$ through the super-exchange process\cite{petta2005coherent,hanson2007spins}. Thus, the ground state becomes a low-spin state for large enough $dV$. We note that the range of $dV$ that we are sweeping is smaller than the Hubbard interactions (on the order of meVs); therefore, the transition is not caused by a direct doublon formation in the detuned site. In addition, the range of detuning over which the high-spin ground state survives is larger than the hybridization $\sim40\,\mu$eV, consistent with the experiment\cite{dehollain2019nagaoka}. This can be reflected by the excited-state spectrum in Fig.~\ref{fig:detuningExcitedState}: the transition between low-spin and high-spin states occurs ``adiabatically'' between the ground states of each section. The Nagaoka gap is always much smaller than the level spacing, which is roughly reflected by the gap between the high-spin ground and excited states.

On the other hand, it is much easier to empty a site compared to doubly occupying one, in a hole-doped system (with three electrons on four sites): the detuning potential only has to compensate the kinetic energy instead of interaction energy to achieve the former. Thus, with a negative $dV$, the $E-dV$ slope flattens out rapidly, except for a small influence from the presence of hybridization. This means that increasing the site energy causes the emptying of the particular dot. For large enough $-dV$, the many-body system becomes an effective empty site plus three singly-occupied dots, or equivalently, a half-filled open-boundary array. Without the ``mobile'' hole in the ``half-filled'' system, the ground state becomes a low-spin state instead of the Nagaoka FM state.

\begin{figure}[!t]
\begin{center}
\includegraphics[width=\columnwidth]{Fig_detuneExctState_D210.jpg}
\caption{\label{fig:detuningExcitedState} (a) The first three excited-state energies in high-spin (blue) and low-spin (red) sectors. The arrow denote the region of Nagaoka phase. (b) The Nagaoka gap and (c) an enlarged energy devolution for the dashed boxed region in (a).
}
\end{center}
\end{figure}

\begin{figure*}[!th]
\begin{center}
\includegraphics[width=2\columnwidth]{Fig_angleTwist_new.jpg}
\caption{\label{fig:angleGroundState} (a) Schematic illustrating the bond rotation in a $2\times2$ system. (b) The ground-state energy for the entire system as a function of the rotating angle $\theta$, calculated for various distances $d$. (c) The Nagaoka gap and (d) the first three excited-state energies in high-spin (blue) and low-spin (red) sectors for $d=210$\,nm. The arrow denotes the region of Nagaoka phase.
}
\end{center}
\end{figure*}

The effect of hybridization is made clear by the dot distance $d$ comparison in Fig.~\ref{fig:detuningGroundState}(b). With increasing distance, the slopes tend towards 0 for negative $dV$ and 1 for positive $dV$, since the increase in distance effectively suppresses any hybridization effects. Interestingly, although the Nagaoka gap decays rapidly for $d>210$\,nm in Fig.~\ref{fig:distDep}, it does not reflect the robustness against potential detuning. In fact, the range of $dV$ where the Nagaoka phase is retained is similar for $d=210$ and 220\,nm. Only after 220\,nm, the range starts to shrink. This is because the robustness of the Nagaoka phase not only depends on the absolute energy gap, but its relative strength compared to the effective tunneling $t$. The fact that $t$ drops by a factor of 2 from $d=210$\,nm to $d=220$\,nm compensates the reduction of the absolute Nagaoka gap.

\subsection{From a Plaquette to a Chain}

By increasing the distance between two of the dots in the plaquette, we can study the four-dot system under different topologies. The Nagaoka theorem applies to a 2D system with periodic boundary conditions. In contrast, a 1D open-boundary system must obey the Lieb-Mattis theorem, which restricts the ground-state solution to the lowest spin sector\cite{lieb1962lieb}. We can gradually change the topology, from a plaquette to a chain, by increasing the angle $\theta$ between two edges in the 2$\times$2 system, as shown in Fig.~\ref{fig:angleGroundState}(a). We again focus on the $d=210$\,nm system first. As shown in Fig.~\ref{fig:angleGroundState}(b), the ground state soon becomes a low-spin state for a rotation angle of $\sim$0.3$^\circ$ [also see Fig.~\ref{fig:angleGroundState}(c)]. The rapid increase of the ground-state energy indicates its sensitivity to the angle, or the topology. This sensitivity can be understood from the excited-state spectrum. The original plaquette has a $C_4$ rotational symmetry, leading to a rotational symmetric ground state. The first and second excited states correspond to the eigenstates of rotation with a factor of $e^{i\pi/2}$ and $e^{-i\pi/2}$, which are degenerate for $\theta=0$. Thus, the extent to which the system ceases to be 2D can be reflected by the energy splitting of these two excited states. As shown in Fig.~\ref{fig:angleGroundState}(d), these two lowest excited states soon separate from each other and the separation becomes comparable with the gap to the ground state for $\theta\sim 0.5^\circ$. This phenomenon indirectly reflects the fact that the system, including its ground state, becomes more like 1D in contrast to 2D, resulting in an $S=1/2$ instead of $S=3/2$ ground state.

Interestingly, the transition from high- to low-spin ground state occurs at very small angles, far before the system becomes 1D geometrically. As Mattis has pointed out, the Lieb-Mattis theorem holds only for a strictly 1D open-boundary system \cite{mattis2003eigenvalues}. That being said, there should be additional mechanisms accounting for the drop of Nagaoka ferromagnetism. The answer to this question might come from the intuition that Nagaoka ferromagnetism is a consequence of constructive interference between the paths that the hole can take through the plaquette, lowering the kinetic energy in the presence of $C_4$ rotational symmetry. This interference effect is quickly lost at even small values of $\theta$, with the broken rotational symmetry leading to unbalanced $x$- and $y$-direction hopping.

Alternatively, one can look at the above reasoning in terms of translational symmetry. Once the hopping between any neighboring sites is dramatically weakened, the system behaves more like an open-boundary chain describable by the Lieb-Mattis theorem. In this sense, the high- to low-spin transition is caused by unbalanced tunneling in the system, rather than geometry. In the experiment, the geometric modification of the system is achieved by tuning the gate potentials, which has a combined effect of increasing the potential barrier between the dots, as well as increasing their separation.

We also examine the transition for different distances $d$, as shown in Fig.~\ref{fig:angleGroundState}(b). Here we observe that the Nagaoka high-spin state is almost equally robust as a function of distance. This can be attributed to the fact that the intrinsic interaction and tunneling scales are almost unchanged when one rotates the two edges, especially for larger distances where the hybridization is negligible. In the former case, only the relative values of the tunneling strengths show up in the path interference, which depends on the rotation angle instead of the absolute tunneling strengths.

\section{Discussion and Conclusion}\label{sec:conclusion}
We described a theoretical, \emph{ab initio} analysis of a quantum dot plaquette system, in which we obtained quantitative agreement with the recent experimental study of the emergence of quantum magnetism through the Nagaoka mechanism. Our work provides theoretical support for the experimentally observed robustness of the Nagaoka state against perturbations such as distance between the dots and potential detuning. Interestingly, one can also find good agreement between experiments and a single-band extended Hubbard model by properly choosing model parameters\cite{dehollain2019nagaoka}. The effective single band should be understood as being comprised of a linear superposition of single-particle electronic orbitals determined by strong inter-orbital interactions. This phenomenological approach, however, has very limited predictive power as it fully relies on fitting parameters to experimental measurements. Our analysis demonstrates that \emph{ab initio} calculations are possible for experimentally relevant systems and can be used to study phenomena beyond the single-band model\cite{nelson2004odd,wu2008orbital,kim2009phase}. Even for the quantitative modeling on a single-band level, we expect the ``bottom-up'' approach to be more accurate than fitting to experimental data. Current experiments can only provide limited information about the excited states and gap sizes, even with the state-of-the-art experimental techniques, and do not allow to determine all parameters of the effective model. With a practical down-folding to the fewer-orbital models, one can further extend the calculation approximately to much larger quantum-dot systems.

Thus, with the focus on a tunable quantum-dot system, we have introduced the \emph{ab initio} exact diagonalization approach, which can be in general applied to different types of artificial quantum simulators. The computational complexity for the model parameter evaluation scales polynomially with the number of sites and orbitals. Calculating the expensive two-center integrals is most costly in the plaquette system. The next level of complexity for these calculations would consider multiple and inhomogeneous Gaussian decompositions, which are significant for stronger hybridized systems or higher-order corrections. These issues have been overcome in modern quantum chemistry using composite atomic basis. Through appropriate fitting using an extended Gaussian basis, we expect to solve these issues by the same means. In any case, the bottleneck of the \emph{ab initio} calculation comes from the bottom-level one-center and two-center integrals Eqs.~\eqref{eq:oneCenterIntegral} and \eqref{eq:twoCenterIntegral}. They have been shown to be efficiently accelerated using GPU-based programming, which can also be directly ported into our systems.

The evaluation of many-body model parameters through the \emph{ab initio} calculation has achieved the goal of precisely modeling an artificial electronic system. Although we here adopt the four-well system and the Nagaoka transition as an example of our approach, motivated by the recent quantum-dot experiment, we would like to emphasize that the \emph{ab initio} exact diagonalization approach can be applied to larger quantum-dot systems with necessary numerical improvements. Unlike the traditional mean-field approaches, a many-body numerical solver like exact diagonalization is always necessary to obtain the ground-state or excited-state wavefunctions. This step is relatively cheap in the current example, but scales up exponentially with the number of sites and electrons. To simulate a larger system, a proper separation of scales might be necessary. For example, if the electron occupation is large, the ``fully occupied'' low energy states may be treated by mean-field theory as a pseudo-potential, to limit the complexity to the bands near the Fermi level. Additionally, the efficiency of the modeling may be further increased employing other many-body numerical approaches including quantum Monte Carlo, density matrix renormalization group, embedding theory, and quantum cluster methods, depending on the purpose of calculation.

Focusing specifically on quantum dot simulators, the accessibility of multiple orbitals and precise treatment of electron interactions could enable a direct simulation of many-body states. Owing to the tunability and measurability of electronic configurations, the quantum dots have been shown to emulate artificial chemical molecules with dominant 2D geometry. For example, the four-dot system investigated in this work can be regarded as an H$_4$ molecule, which is a standard platform for testing correlated quantum chemistry methods. Hence the quantum-dot simulators can be used to find the many-body electron state in a Born-Oppenheimer assumption.

Looking beyond quantum dot systems, this approach can be naturally extended to Rydberg atoms or cold molecules by replacing the Coulomb interaction $W(\rbf_1-\rbf_2)$ with the Lennard-Jones potential and making $V(\rbf)$ a standing-wave potential. The breaking of rotational symmetry in $V(\rbf)$ may cause more computational complexity, which can be overcome using some of the efficient integration implementations mentioned above. Moreover, the majority of the optical lattice studies concern bosons. The \emph{ab initio} exact diagonalization framework can be applied to bosonic systems by replacing the fermionic basis states represented by Slater determinants with bosonic product states represented by permanents. In general, this approach holds the promise to push the boundaries of predictability and quantitative accuracy in the ever-expanding zoo of quantum simulators that are being implemented.

\section*{Acknowledgements}
We thank for insightful discussions with S. Das Sarma and B. Wunsch. This work was supported at Harvard University by the Harvard-MIT Center for Ultracold Atoms, NSF Grant No.~DMR-1308435, the AFOSR-MURI Photonic Quantum Matter (Award No.~FA95501610323), and the DARPA DRINQS program (Award No.~D18AC00014). Y.W. was supported by the Postdoctoral Fellowship in Quantum Science of the Harvard-MPQ Center for Quantum Optics. L.M.K.V. thanks the NSF-funded MIT-Harvard Center for Ultracold Atoms for its hospitality. J.P.D., U.M., and L.M.K.V. acknowledge grants from the Netherlands Organisation for Scientific Research (FOM projectruimte and NWO Vici). M.R. thanks the Villum Foundation for support. This research used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility operated under Contract No. DE-AC02-05CH11231.