Ab Initio Exact Diagonalization Simulation of the Nagaoka Transition in Quantum Dots

TL;DR

This paper introduces an ab initio exact diagonalization method to model correlated electrons in quantum dots, accurately studying the Nagaoka transition and ferromagnetism, aligning well with experimental data and advancing quantum simulator analysis.

Contribution

The work presents a novel ab initio exact diagonalization framework that includes realistic multi-orbital effects to analyze ferromagnetism and the Nagaoka transition in quantum dots.

Findings

Good agreement with experimental measurements of magnetism.

Robustness of ferromagnetic states under various conditions.

Insights into multi-orbital effects on the Nagaoka transition.

Abstract

Recent progress of quantum simulators provides insight into the fundamental problems of strongly correlated systems. To adequately assess the accuracy of these simulators, the precise modeling of the many-body physics, with accurate model parameters, is crucially important. In this paper, we introduce an \emph{ab intio} exact diagonalization framework to compute the correlated physics of a few electrons in artificial potentials. We apply this approach to a quantum-dot system and study the magnetism of the correlated electrons, obtaining good agreement with recent experimental measurements. Through dot separation, potential detuning and control of single tunneling, we examine the Nagaoka transition and determine the robustness of the ferromagnetic state. While the standard Nagaoka theorem considers only a single-band Hubbard model, in this work we perform extensive $ab$ $intio$ calculations that include realistic multi-orbital conditions in which the level splitting is smaller than the interactions. This simulation complements the experiments and provides insight into the formation of ferromagnetism in correlated systems. More generally, our calculation sets the stage for further theoretical analysis of analog quantum simulators at a quantitative level.

Tables

Table 1. Table 1: Effect of system parameters (definition of these parameters can be found in the Sec.  III.2 ) on ground-state energies and the Nagaoka gap obtained by various models for d = 210 𝑑 210 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.

Model	Ground-State Energy	Nagaoka Gap
$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

Equations

$$\left[-\frac{{\hslash }^2}{2m_e^{\star }}\left(\frac{{\partial }^2}{\partial r^2}+\frac{1}{r}\frac{\partial }{\partial r}+\frac{1}{r^2}\frac{{\partial }^2}{\partial {\varphi }^2}\right)+V(r)\right]\psi (\text{r})=E\psi (\text{r}),$$

$$\psi (\text{r})=\frac{\chi (r)}{\sqrt{2\pi r}}{e}^{im\varphi }=\frac{\chi (r)}{\sqrt{r}}\phi (\varphi ),$$

$${\int }_0^{\infty }{\chi }_n(r{)}^{\ast }{\chi }_{n^{\prime }}(r)dr={\delta }_{n{n}^{\prime }}.$$

$$-\frac{{\hslash }^2}{2m_e^{\star }}\left[\frac{d^2\chi }{d{r}^2}-\frac{m^2-1/4}{r^2}\chi \right]+V(r)\chi =E\chi .$$

$$\psi_{nm}(\textbf{r})=\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.$$

$$\mathcal{H}=\sum _a^{N_e}\left[-\frac{{\hslash }^2}{2m_e^{\star }}{\nabla }_a^2+\sum _i^{N_{\mathrm{well}}}V({\text{r}}_a-{\text{R}}_i)\right]+\sum _{a\ne b}\frac{e^2}{4\pi ϵ|{\text{r}}_a-{\text{r}}_b|},$$

$${\tilde{\psi }}_{\mu \sigma }(\text{r})=\sum _{\nu }{X}_{\nu \mu }{\psi }_{\nu \sigma }(\text{r}).$$

$$\int d{\text{r}}^3{\psi }_{\mu \sigma }^{\ast }(\text{r}){\psi }_{\nu {\sigma }^{\prime }}(\text{r})=S_{\mu \nu }{\delta }_{\sigma {\sigma }^{\prime }}.$$

$$\int d{\text{r}}^3{\tilde{\psi }}_{{\mu }_1{\sigma }_1}^{\ast }(\text{r}){\tilde{\psi }}_{{\mu }_2{\sigma }_2}(\text{r})=X^{†}SX{\delta }_{{\sigma }_1{\sigma }_2}.$$

$$|\{i_k{\alpha }_k{\sigma }_k\}⟩=c_{i_N{\alpha }_N{\sigma }_N}^{†}\cdots c_{i_2{\alpha }_2{\sigma }_2}^{†}{c}_{i_1{\alpha }_1{\sigma }_1}^{†}|0⟩,$$

$$\mathcal{H}={\mathcal{H}}_{\mathrm{non}-\mathrm{int}}+{\mathcal{H}}_{\mathrm{int}}.$$

$${\mathcal{H}}_{\mathrm{non}-\mathrm{int}}=\sum _{i\alpha \sigma }\sum _{j\beta \sigma }{t}_{\genfrac{}{}{0ex}{}{ij}{\alpha \beta }}{c}_{i\alpha \sigma }^{†}{c}_{j\beta \sigma }.$$

$$h_{\mu \nu }$$

$$t_{\mu \nu }$$

$${\mathcal{H}}_{\mathrm{int}}=\sum _i{\mathcal{H}}_i^{(\mathrm{OS})}+\sum _{ij}{\mathcal{H}}_{ij}^{(\mathrm{LR})}.$$

$${\mathcal{H}}_i^{(\mathrm{OS})}$$

$${\Xi }_{\genfrac{}{}{0ex}{}{\mu {}_1\nu {}_1}{\mu {}_2\nu {}_2}}=\iint d{\text{r}}_1^{d}d{\text{r}}_2^{d}W(|{\text{r}}_1-{\text{r}}_2|){\psi }_{{\mu }_1}({\text{r}}_1{)}^{\ast }{\psi }_{{\mu }_2}({\text{r}}_2{)}^{\ast }{\psi }_{{\nu }_1}({\text{r}}_1){\psi }_{{\nu }_2}({\text{r}}_2).$$

$$U_{{\alpha }_1{\alpha }_2}^{\prime }$$

$$J_{{\alpha }_1{\alpha }_2}$$

$${\mathcal{H}}_{ij}^{(\mathrm{LR})}$$

$$V_{\alpha \beta }^{ij}$$

$$K_{\alpha \beta }^{ij}$$

$$V_{\alpha \beta }^{ij\prime }$$

$$K_{\alpha \beta }^{ij\prime }$$

$${\mathcal{H}}_{\mathrm{int}}=\frac{1}{2}\sum _{i_1{\alpha }_1{\sigma }_1}\sum _{i_2{\alpha }_2{\sigma }_2}\sum _{j_1{\beta }_1{\sigma }_1^{\prime }}\sum _{j_2{\beta }_2{\sigma }_2^{\prime }}\mathcal{W}(j_1,{\beta }_1,{\sigma }_1^{\prime };j_2,{\beta }_2,{\sigma }_2^{\prime }|i_2,{\alpha }_2,{\sigma }_2;i_1,{\alpha }_1,{\sigma }_1)c_{j_1{\beta }_1{\sigma }_1^{\prime }}^{†}{c}_{j_2{\beta }_2{\sigma }_2^{\prime }}^{†}{c}_{i_2{\alpha }_2{\sigma }_2}{c}_{i_1{\alpha }_1{\sigma }_1}.$$

$$\mathcal{W}(j_1,{\beta }_1,{\sigma }_1^{\prime };j_2,{\beta }_2,{\sigma }_2^{\prime }|i_2,{\alpha }_2,{\sigma }_2;i_1,{\alpha }_1,{\sigma }_1)$$

$${\mathcal{H}}_{\mathrm{int}}$$

$$U_{ij}^{\sigma {\sigma }^{\prime }}({\beta }_1,{\beta }_2|{\alpha }_2,{\alpha }_1)=U_{ji}^{{\sigma }^{\prime }\sigma }({\beta }_2,{\beta }_1|{\alpha }_1,{\alpha }_2)$$

$$J_{ij}^{\sigma {\sigma }^{\prime }}({\beta }_1,{\beta }_2|{\alpha }_2,{\alpha }_1)=J_{ji}^{{\sigma }^{\prime }\sigma }({\beta }_2,{\beta }_1|{\alpha }_1,{\alpha }_2).$$

$${\mathcal{H}}_i^{(\mathrm{OS})}$$

$$U_{\alpha }$$

$$=\iint d{\text{r}}_1^{d}d{\text{r}}_2^{d}W(|{\text{r}}_1-{\text{r}}_2|)|{\tilde{\psi }}_{i\alpha \sigma }({\text{r}}_1){|}^2|{\tilde{\psi }}_{i\alpha \bar{\sigma }}({\text{r}}_2){|}^2.$$

$$U_{{\alpha }_1{\alpha }_2}^{\prime }-J_{{\alpha }_1{\alpha }_2}$$

Full text

I Introduction

Strong correlations are at the heart of many important phenomena in condensed matter systems, including unconventional superconductivityBednorz and Müller (1986), quantum magnetismFazekas (1999) and fractional quantum Hall statesStormer and Tsui (1983). These phenomena have a wide range of applications in material design, energy science, and quantum informationDagotto (1994). 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 diagonalizationLin (1990), quantum Monte CarloFoulkes et al. (2001), density-matrix renormalization groupSchollwöck (2005), 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 systemsBloch et al. (2008); Lewenstein et al. (2012); Gross and Bloch (2017) and have recently begun exploring fermionic systemsParsons et al. (2015); Greif et al. (2016); Parsons et al. (2016); Boll et al. (2016); Brown et al. (2019); Nichols et al. (2019). 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 solidsVan der Wiel et al. (2002); Hanson et al. (2007); Nielsen and Bhatt (2007); Oguri et al. (2007); Thalineau et al. (2012); Seo et al. (2013); Hensgens et al. (2017); Mukhopadhyay et al. (2018). 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 parametersHilker et al. (2017); Mazurenko et al. (2017); Barthelemy and Vandersypen (2013). 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 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, ab initio calculations have made great progress towards describing systems with increasing complexity. With the help of the Gaussian basisBinkley et al. (1980); Ditchfield et al. (1971); Frisch et al. (1984); Hariharan and Pople (1973); Hehre et al. (1969), the computational cost has been largely reduced, making the simulation of large molecules possible. In addition to the basic Hartree-Fock methodHartree (1928); Hartree and Hartree (1935); Fock (1930a, b); Slater (1930), many advanced post-Hartree-Fock wavefunction-based methods (coupled clusterShavitt and Bartlett (2009), configuration interactionsSherrill and Schaefer III (1999); Pople and Head-Gordon (1987), etc.) and multi-reference methodsHegarty and Robb (1979); Eade and Robb (1981); Yamamoto et al. (1996); Andersson et al. (1990, 1992) 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 systemsUfimtsev and Martinez (2008a, b); DePrince III and Hammond (2011); Hohenstein et al. (2015); Song and Martínez (2018). 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 diagonalizationMaksym and Chakraborty (1990); Merkt et al. (1991); Pfannkuche et al. (1993); Hawrylak and Pfannkuche (1993); Yang et al. (1993); Hawrylak (1993); Palacios et al. (1994); Wojs and Hawrylak (1995); Oaknin et al. (1995); Wojs and Hawrylak (1996); Maksym (1996); Hawrylak et al. (1996); Wojs and Hawrylak (1997a, b); Eto (1997); Maksym (1998); Reimann et al. (2000); Mikhailov (2002) and quantum Monte CarloBolton (1996); Harju et al. (1999a, b, 2002); Siljamäki et al. (2002). The simplified treatment was also extended to double-dot systemsBurkard et al. (1999); Li et al. (2009); Nielsen et al. (2012); Stepanenko et al. (2012); Yannouleas et al. (2016); Brandt et al. (2017). 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 statesAbolfath and Hawrylak (2006); Stopa et al. (2006); Abolfath et al. (2006). However, these are the largest quantum-dot systems that have been subject to 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 modelsStafford and Sarma (1994, 1997); Stafford et al. (1998); Kotlyar et al. (1998a, b).

A recent experiment Dehollain et al. (2019) 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 ferromagnetismNagaoka (1966). The success of the experiment in Ref. Dehollain et al., 2019 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 et al Dehollain et al. (2019). 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 ab initio inputs is also necessary.

For both of these purposes, we hereby introduce an 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 experimentsDehollain et al. (2019). 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 experimentDehollain et al. (2019).

The description of this model and calculation will gradually increase in complexity. In Sec. II, 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. III. By adjusting the model to represent a four-well system, in Sec. IV we then explore the quantum magnetism and especially the Nagaoka transition using the ab initio exact diagonalization approach. Finally, we conclude and discuss the future directions of our approach in Sec. V.

II Single Electron in a Single-Quantum Well

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 dotpot . 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(\textbf{r})=-V_{0}e^{-{|\textbf{r}|^{2}}/{2\sigma}}$ in this paper [see Fig. 1(a)], as a typical description of the finite-size quantum dotAbolfath et al. (2006); Abolfath and Hawrylak (2006). Here $\textbf{r}=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 Schrödinger equation is

[TABLE]

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

[TABLE]

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

[TABLE]

Then we obtain the radial differential equation

[TABLE]

It can be numerically solved using the finite difference approximation. Choosing the angular part being real for numerical convenience, we define the single-well wavefunction as

[TABLE]

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. 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 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$ (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. 1, resulting in the level splitting of the degenerate statesvan Beveren et al. (2005). However, as we will show in Sec. IV, 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.

III Many-Body Model

With multiple potential wells, the general Hamiltonian for a many-body system among $N_{\rm well}$ wells is

[TABLE]

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(\textbf{r})$ 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 theorySzabo and Ostlund (2012), we construct the basis using a superposition of the single-well wavefunctions

[TABLE]

For simplicity in notation, we collapse the coordinate and orbital indices as $\mu=(i,\alpha)$, and denote $\psi_{\mu\sigma}(\textbf{r})=\psi_{\alpha}(\textbf{r}-\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 Schrödinger equation. The overlap matrix among different single-well wavefunction basis is

[TABLE]

Thus, the overlap matrix among the new basis functions is

[TABLE]

By setting the requirement $X^{\dagger}SX=I$ and considering $S$ being positive-definitive, a standard choiceSzabo and Ostlund (2012) is $X=S^{-1/2}$. This selection results in a new orthonormal basis set $\{\tilde{\psi}_{\mu\sigma}(\textbf{r})\}$.

Representing the many-body wavefunction in the Fock space, spanned by the Slater determinants of $\{\tilde{\psi}_{\mu\sigma}(\textbf{r})\}$, we obtain the the second quantization of the many-body statesHelgaker et al. (2014):

[TABLE]

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

[TABLE]

The non-interacting part ${\mathcal{H}}_{\rm non-int}$ corresponds to the hopping of an electron across orbitals and wells, while the interacting part ${\mathcal{H}}_{\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.

III.1 Non-Interacting Part of Hamiltonian

The non-interacting part of the Hamiltonian is quadratic in fermionic operators

[TABLE]

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

[TABLE]

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. (7) into (13), resulting in

[TABLE]

Here $t_{\mu\nu}$ defines the site energy (diagonal) and hybridization (off-diagonal) of the orthonormal orbitals, which appears in Eq. (11). 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.

III.2 Interacting Part of Hamiltonian

As a typical choice in condensed matter, we restrict the interaction part ${\mathcal{H}}_{\rm int}$ to four-fermion termsSzabo and Ostlund (2012). As specified in Eq. (21), 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. (22). 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 processesres , 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. 2]. 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 levelsnon . 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 stateDagotto (2013).

After these simplifications, the interacting part of the Hamiltonian can be decomposed as the on-site and (two-center) long-range parts

[TABLE]

The standard derivations of each term in the interacting Hamiltonian is present in Appendix A. The on-site interaction Hamiltonian can be written as

[TABLE]

This is the known as the multiplet model, widely used to describe the valence electrons in the transition metal systemsDagotto (2013). The corresponding scattering processes are sketched in Fig. 3. For convenience, the site index is removed on the right-hand side, while in an inhomogeneous system [such as the modulations in Sec. IV], one should consider it specifically for each individual site. Due to a symmetry consideration [see discussions in Appendix B], it is usually convenient to calculate the interaction parameters using the original single-well basis obtained from Eq. (5), through

[TABLE]

Here $W(|\textbf{r}_{1}\!-\!\textbf{r}_{2}|)={e^{2}}/{4\pi\epsilon|\textbf{r}_{1}-\textbf{r}_{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. (7), we have

[TABLE]

These parameters define the on-site interactions among the orthonormal orbitals.

Similarly, the long-range interactions are written as

[TABLE]

As sketched in Fig. 4, $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

[TABLE]

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.  IV), 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.

IV Simulation of Four-Well Quantum Dot System: Probing Nagaoka Magnetism

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 NagaokaNagaoka (1966), 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. 5, 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 limitvon Stecher et al. (2010).

A recent experiment has studied Nagaoka magnetism using a quantum dot array in a $2\times 2$ plaquette configurationDehollain et al. (2019). 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 ab initio exact diagonalization approach described above to extract the precise many-body model of the $2\times 2$ quantum dot plaquette and quantitatively reproduce the Nagaoka conditions that were explored with the experimental system.

IV.1 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 nmMukhopadhyay et al. (2018); Dehollain et al. (2019). 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}^{2}m_{e}^{\star}}\approx 0.114$ meV. Applying this scale to the eigen-spectrum solved in Fig. 1 (i.e., $V_{0}=100\,{\hbar^{2}}/{a_{0}^{2}m_{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.

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}\approx 2.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 6 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. (16). As shown in Table 1, only the long-range Coulomb interaction $V$ obviously affects the ground-state energy, by order of 1 meV, while others contribute to $\sim 0.01$ meV. However, as stated before and now made clear in Table 1, 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. (22). 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. Dehollain et al., 2019 is indeed caused by a Nagaoka-like mechanism.

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 ($\sim 0.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 problemcoh . 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 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 7 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$eVMukhopadhyay et al. (2018); Dehollain et al. (2019), which in the 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 ab initio calculation using only very limited experimental input match quantitatively with the experiment in Ref. Dehollain et al., 2019. Therefore, we believe the 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 ignoredcon . We perform exact diagonalization to solve this 60-orbital spinful system, using the parallel Arnoldi approachLehoucq et al. (1998); Jia et al. (2018).

IV.2 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. 8(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 HubbardNagaoka (1966); Mattis (2003) and extended Hubbard modelsKollar et al. (1996). Here we show its validity in a multi-orbital system.

As shown in Fig. 8(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. Dehollain et al., 2019 through a comparison between experimentally measured parameters with a fitted single-band model.

IV.3 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. 9(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 ab initio calculation.

The results from this study, shown in Fig. 9(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 experimentDehollain et al. (2019). 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 processPetta et al. (2005); Hanson et al. (2007). 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 $\sim 40\,\mu$eV, consistent with the experimentDehollain et al. (2019). This can be reflected by the excited-state spectrum in Fig. 10: 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.

The effect of hybridization is made clear by the dot distance $d$ comparison in Fig. 9(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. 8, 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.

IV.4 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 sectorLieb (1962). 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. 11(a). We again focus on the $d=210$ nm system first. As shown in Fig. 11(b), the ground state soon becomes a low-spin state for a rotation angle of $\sim$0.3∘ [also see Fig. 11(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. 11(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 Mattis (2003). 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. 11(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.

V Discussion and Conclusion

We described a theoretical, 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 parametersDehollain et al. (2019). 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 ab initio calculations are possible for experimentally relevant systems and can be used to study phenomena beyond the single-band modelNelson et al. (2004); Wu (2008); Kim et al. (2009). 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 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 ab initio calculation comes from the bottom-level one-center and two-center integrals Eqs. (B) and (35). 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 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 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 H4 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(\textbf{r}_{1}-\textbf{r}_{2})$ with the Lennard-Jones potential and making $V(\textbf{r})$ a standing-wave potential. The breaking of rotational symmetry in $V(\textbf{r})$ 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 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.

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.

Appendix A Derivation of the Interacting Part of the Hamiltonian

Restricting to the four-fermion terms, the second-quantized Hamiltonian can be generically expressed asSheng et al. (2005)

[TABLE]

Substituting the wavefunctions into it, we obtain

[TABLE]

Using the simplification mentioned in the main text, the interaction terms can be categorized into $\sigma_{1}=\sigma_{1}^{\prime}$ and $\sigma_{1}=\sigma_{2}^{\prime}$ parts. Denoting these two parts as $U$ and $J$, we obtain

[TABLE]

Note the $U$ and $J$ terms are not completely independent, since $U_{ij}^{\sigma\sigma}\equiv J_{ij}^{\sigma\sigma}$. Additionally, we also have the permutation symmetry

[TABLE]

For electrons on a single lattice site, the generic form in Eq. (23) reduces to

[TABLE]

Note that the spin-parallel Hund term $J_{\alpha_{1}\alpha_{2}}$ is the same as the spin-parallel Hubbard term $U^{\prime}_{\alpha_{1}\alpha_{2}}$ with a sign flip.

Therefore, the on-site Hubbard interaction in Eq. (A) is

[TABLE]

The Hubbard interaction is dominant among the interaction terms due to the maximal overlap of wavefunctions. The remaining terms in a single-well interaction are all the inter-orbital interactions. The spin-parallel interaction is

[TABLE]

while spin-anti-parallel interaction is

[TABLE]

Given that the two-body interaction $W(|\textbf{r}_{1}-\textbf{r}_{2}|)$ (typically Coulomb) does not involve spin degrees of freedom, the first term of Eq. (27) is equal to the anti-parallel spin contribution in Eq. (28). Naturally, one can split the entire parallel spin interactions in Eq. (27) into charge and Hund’s part by assuming $U^{\prime}_{\alpha_{1}\alpha_{2}}=\bar{U}^{\prime}_{\alpha_{1}\alpha_{2}}$. This partition also guarantees the equivalence of the two exchange coefficients

[TABLE]

Therefore, we obtain the on-site interacting Hamiltonian Eq. (16) in the main text.

Then following Eq. (19) in the main text, we can evaluate the interaction parameters in the long-range part of Hamiltonian ${\mathcal{H}}_{ij}^{\rm(LR)}$. Similar to the on-site terms, the “off-diagonal” terms of $V_{\alpha\beta}$ and $V^{\prime}_{\alpha\beta}$ are absorbed by the corresponding exchange terms for parallel spins. Therefore, we can write expressions for each of the relevant long-range terms

[TABLE]

and transform them to the original basis, resulting in Eq. (III.2) in the main text.

Appendix B Evaluation of the Integrals

The algebraic representations of the model parameters Eqs. (18) and (III.2) concentrate all integration calculations in the evaluation of $\Xi_{\mu_{1}\nu_{1}\atop\mu_{2}\nu_{2}}$ in the single-well basis. This evaluation is not trivial, since the direct expression Eq. (17) contains a $2\times d$-dimensional integral with singularities, which cannot be computed efficiently even with supercomputersHelgaker et al. (2014). However, taking advantage of the rotational symmetry of the quantum well, the calculation can be significantly simplified.

Let us first look at the dominant part – the one-center integral, where all four wavefunctions are centered in the same well. Taking advantage of the rotational invariance of $W(|\textbf{r}_{1}\!-\!\textbf{r}_{2}|)$, one can simplify the integral through the Wigner-Eckart theorem. Specifically, for a Coulomb-type interaction, we have the Laplacian expansion

[TABLE]

in which $P_{l}(x)$ is the Legendre polynomial. With $W(|\textbf{r}_{1}\!-\!\textbf{r}_{2}|)=e^{2}/4\pi\epsilon|\textbf{r}_{1}\!-\!\textbf{r}_{2}|$, we can simplify the one-center integral to a sequence of two-dimensional integrals as

[TABLE]

then Eq. (17) is expanded as

[TABLE]

Note, that $\Theta_{\mu_{1}\nu_{1}\atop\mu_{2}\nu_{2}}^{(l)}$ is symmetric under exchange of 1 and 2 indices, while $\Xi_{\mu_{1}\nu_{1}\atop\mu_{2}\nu_{2}}$ is usually not symmetric except in special cases where $\{\mu_{1},\nu_{1}\}=\{\mu_{2},\nu_{2}\}$. The integral decays rapidly with the increase of $l$. With fine enough spatial grids and angular momentum truncation, the one-center integral can be evaluated up to machine precision.

In contrast, the two-center integral involves more computational complexity. Here, rotational symmetry is not maintained; therefore, there is no direct separation of variables. However, we know that the ground state and the norm of low-lying excited-state wavefunctions can be well estimated by different Gaussian functions. This provides a way to estimate the density-density correlation among the two-center integrals. If the density distribution is written as

[TABLE]

the two-center integral can be decomposed in the center-of-mass frame

[TABLE]

Now, the integral is reduced to a two-dimensional integral in the reduced coordinates $\Delta\textbf{r}$, which can be solved by using the Riemann integral or the Laplacian expansion as mentioned above. Note that the Gaussian integral provides only an estimation of the realistic two-center interaction. A more precise treatment involves the decomposition of multiple Gaussian bases and its derivativesMcMurchie and Davidson (1978), which forms the foundation of electronic structure theory and is beyond the scope of this work.