Theory of charge stability diagrams in coupled quantum dot qubits
Nathan L. Foulk Sankar Das Sarma


    
    
    
    Condensed Matter Theory Center and Joint Quantum Institute, Department of Physics, University of Maryland, College Park, Maryland 20742-4111 USA
§ ABSTRACT 
 We predict large regions of the charge stability diagram using a multi-band and multi-electron configuration interaction model of a double quantum dot system.
 We account for many-body interactions within each quantum dot using full configuration interaction and solve for single-particle density operators.
 This allows charge states to be predicted more accurately than the extensively used classical capacitance model or the single-band Hubbard model. 
 The resulting single-particle mixed states then serve as inputs into an atomic orbital picture that allows for the explicit calculation of the underlying Hubbard model parameters by performing the appropriate integrals. This numerical approach allows for arbitrary choices of electrostatic potential and gate geometry.
 A common assumption when calculating charge stability diagrams from the Hubbard model is that the charge stability diagrams are periodic, but we find that the tunnel couplings for valence electrons in dots with N=3 electrons are significantly enhanced when compared to single-electron dots. 
 This difference is apparent in the charge stability diagram for higher occupancy Coulomb diamonds. 
 We also quantitatively explore how the barrier gate strength and dot pitch impact this behavior.
 Our work should help improve the future realistic modeling of semiconductor-dot-based quantum circuits.

Theory of charge stability diagrams in coupled quantum dot qubits
    S. Das Sarma
    September 9, 2024 – Version 1.0
=================================================================



§ INTRODUCTION

Quantum dot based semiconductor spin qubits <cit.> are a promising candidate for scalable quantum computing, benefitting from long coherence times <cit.>, fast gate operations  <cit.>, and integration possibilities with the existing semiconductor industry <cit.>.
Using lithographically defined metal gates, quantum dots are formed in the 2D quantum well of a semiconductor heterostructure.
These dots are extremely small–on the order of 100 nm–and their small size provides much of the promise as one of the few qubit implementations that are serious candidates for a fully error-corrected and scalable quantum computer <cit.>—the scalability because of the compatibility with the vast semiconductor electronics industry.
Each dot is loaded with a predetermined number of electrons, and the qubits are represented by the spin states of these electrons.
Single spin manipulation can be performed with electrical or magnetic means and the Heisenberg exchange interaction allows for fast and predictable electrical two-spin control between nearest neighbors <cit.>.

The simplest qubit encoding is the Loss-Divincenzo qubit, where the two-level system is a single electron spin confined to a single quantum dot <cit.>.
Many other qubit encodings exist, including singlet-triplet qubits <cit.>, exchange-only qubits <cit.>, and the flopping-mode qubit <cit.> to name a few.
However, all encoding schemes require precise control over the electron number of each dot (although a restriction to one electron per dot is not essential for qubit operations) <cit.>.
Therefore, one must tune each device so that each dot contains the correct electron number, and this must be done each time the device is cooled and initialized.
In addition, this number of electrons must be preserved during the gate operations—a 1-electron dot cannot suddenly become a 3- or 5-electron dot due to local fluctuations.
A typical spin qubit device contains multiple dots arranged in a linear array using an alternating pattern of plunger gates and barrier gates.
An example of such a device is found in Fig. <ref>(b).
Each quantum dot is electrostatically defined under a plunger gate.
Plunger gate voltages are raised or lowered to control the electron number of each dot.
The two dots are separated by a distance known as the pitch.
Barrier gates (between the dots) are primarily used to modulate the wavefunction overlap between two dots.
The gate voltages on the plunger and barrier gates are the most important variables in determining the charge state of the dots.
Modern quantum dot devices are made of many more gates and components than just these two, but the plunger and barrier gates occupy a central role.
Any minimal model must include the role of these two gates, one controlling each dot and the other controlling the inter-dot coupling.

Charge stability diagrams are used to visualize and map out the different charge regimes in a double quantum dot, thus enabling the precise control over electron occupancy necessary for qubit operations.
The plunger gate voltages for two neighboring dots are varied, and the conductance is measured. 
A nonzero conductance signifies that electrons are tunneling into the device or between dots. 
Modern devices also employ nearby quantum dots for charge-sensing measurements, where the charge-sensor signal measures the charge states capacitively.
An example of a charge stability diagram is shown in Fig. <ref>a.
The diagram typically consists of sharp peaks of transitions and large plateaus of charge state stability.
These different plateaus demarcate different charge states for a double dot system. Each plateau, or cell—referred to as a “Coulomb diamond.”—indicates fixed electron numbers.




These charge stability diagrams are essential for determining not only the initial gate voltages necessary for reliable charge initialization but also for mapping paths across the stability diagram, which are needed to carry out precise gate operations on the qubits.
Recently, charge stability diagrams have also been used as inputs to machine learning approaches for automated control of quantum dot arrays <cit.>.
The direct calculation of realistic charge stability diagrams given the underlying microscopics of a device is therefore a timely need in the quantum computing research community in order to cost-effectively produce qubit designs and to generate input data for machine learning algorithms for autonomous control.

As a complement to direct measurement, one can predict the charge boundaries of the stability diagram. 
This can be done on varying levels of sophistication, the simplest of which is the classical capacitance model of the double dot system <cit.>.
The classical capacitance model uses the physical dimensions of the dots and their pitch to estimate the electrostatic parameters of the system, such as the charging energy of each dot and the Coulomb interaction between the dots.
This approach correctly identifies the charge regimes in the absence of quantum fluctuations. 
In order to include the effects of quantum fluctuations—which are obviously significant for double quantum dot systems—one can use the Hubbard model with nonzero tunneling <cit.>.
The major Hubbard parameters that affect the stability diagram are the onsite Coulomb repulsion U, the interdot Coulomb repulsion V, and the tunnel coupling t.
The generalized Hubbard model reproduces the capacitance model exactly when t = 0. 
When t > 0, however, the sharp corners of the diamonds become rounded, with larger tunnel constants leading to more rounding (see Fig. <ref>c).
The principal difficulty remains in choosing realistic Hubbard parameters for a given physical system. 
Thus, the generalized Hubbard model, introduced in this context in Refs. <cit.>, considerably simplifies charge stability modeling provided that the microscopic parameters t, U, and V are known for the quantum dot system.

In the simplest case, given model confinement potentials and single electron occupancy, the estimation of these effective Hubbard parameters reduces to integrals of the single particle eigenstates <cit.>.
Such an approach is also implicitly based on charge stability diagram periodicity.
One chooses a single point (V_L,V_R) within the (N_L,N_R) = (1,1) cell, calculates the wavefunctions of each dot separately, and estimates the corresponding Hubbard parameters, consistent with a Hund-Mulliken atomic orbital picture.
Those parameters are then extended to the entire stability diagram in a way that produces perfectly periodic Coulomb diamonds.
Experimental stability diagrams superficially justify this approach, as their Coulomb diamonds are often close to periodic with regard to cell size and location. 
However, the smoothness of the cell corners often varies drastically across the diagram <cit.>.
In order to improve this approach of charge stability diagram simulation, calculations featuring multiple electrons in each dot must be included so as to go beyond the (N_L,N_R) = (1,1) Coulomb diamond constraint.
We would expect the Hubbard parameters to depend on device-specific details, such as gate geometry, barrier gate level, electron occupancy, and both magnetic and electrostatic disorder. 
Incorporating such device-specific details requires the use of extensive, system-dependent numerics, rather than analytic approaches.

Going beyond the (N_L,N_R) = (1,1) Coulomb diamond also addresses the practical needs of the spin qubit community, since a common practice is to load a quantum dot with three electrons rather than one (in principle, any odd number of electrons per dot suffices, but three seems to be the optimal number).
The spin state of the valence electron is then used as the qubit's two-level system <cit.>.
Calculating the Hubbard parameters of multi-electron dots is significantly more complex than single-electron dots. 
The presence of other electrons perturbs the valence electron's wavefunction from the simple analytic expressions often used when assuming simplified potentials.
This increased electron number also complicates numerical treatments of the system, for the simple fact that simulating six-electron systems is very computationally expensive, as the relevant Hilbert space grows exponentially with the number of electrons.
An approach that leverages the computational efficiency of the atomic orbital projection method of previous work toward understanding N=3 quantum dots is therefore needed. 
This is our goal in the current work utilizing the configuration interaction method, which is the standard technique in quantum chemistry for calculating the electronic structure of multielectron atoms and molecules.

We first present our approach of using full configuration interaction (FCI) to project an accurate representation of N=3 quantum dots in the effective Hubbard model space for simulating charge stability diagrams.
Separating the double dot system into two effective single dot potentials, we solve for the three electron wavefunction of each dot.
Using the natural orbitals of the one-electron reduced density matrix (1RDM) of each dot, we approximate the atomic orbitals and calculate the relevant Hubbard parameters.
We then present our results, emphasizing the novelties introduced by the many-body treatment, including enhanced valence tunnel coupling and an understanding of how barrier gate strength and device layout impact the system's projection onto the Hubbard model.
Finally, we discuss the implications of our findings for simulating charge stability diagrams and the quantum dot semiconductor spin qubits field at large.
Our explicit results assume the qubit platform to be based on quantum dots on the Si [100] surface, but the theory and the general method are applicable to all semiconductor quantum dot structures.



§ MODEL


We calculate charge stability diagrams of double quantum dot systems using a generalized Hubbard model.


    Ĥ =  μ_1 n_1 +  μ_2 n_2 + ∑_iU/2 n_i (n_i - 1) + V n_1 n_2  
     
     + ∑_m (t_m c^†_1m c^_2m + H.c.)

where n_i is the total electron dot occupancy over all levels for dot i, and c^†_ij is the creation operator for the jth level for the ith dot. 
Many terms can be included in a generalized Hubbard model, but the terms that dominate the charge stability diagram are t, U, and V, which are necessary in the minimal effective model (and should suffice in most situations) since interaction terms beyond nearest-neighbors are likely to be vanishingly small.
All other terms can be completely ignored, and the stability diagram will be almost exactly the same.
This Hamiltonian can be rewritten as 

    Ĥ = ∑_i,j F_ij c_i^† c_j + ∑_i,j,k,l G_ijkl c_i^† c^†_j c_k c_l.

Given the appropriate single particle wavefunctions that correspond to c_i^†|0⟩ and c_j^†|0⟩, where |0⟩ is the vacuum state, we can calculate F_ij

    F_ij = ∫Ψ_i^*(𝐫) ĤΨ_j(𝐫) d 𝐫.

This can easily be extended to two-body terms as well: 

    G_ijkl = ∬Ψ_i^*(𝐫)Ψ_j^*(𝐫') ĤΨ_k(𝐫) Ψ_l(𝐫') d𝐫 d𝐫'.

However, the difficulty lies in extracting these mutually orthogonal, single-particle wavefunctions.
In reality, these single-particle states are entangled, and care must be taken to preserve their orthogonality.
We first calculate the many-body wavefunction.

We utilize full configuration interaction (FCI) <cit.> to calculate the multi-electron ground state. 
These FCI calculations are performed within the effective mass approximation using the following Hamiltonian: 

    H(𝐫) = ∑_i ((𝐩_i^2 + e𝐀(𝐫_i))^2/2m_e^*  + V(𝐫_i) + μ_B 𝐒_𝐢·𝐁)
    
     + ∑_i,je^2/4πϵ1/|𝐫_i - 𝐫_j| .

where m^*_e is the effective mass for silicon in the in-plane directions, approximately m^*_e ≈ 0.19m_0, where m_0 is the electron's rest mass. 
Additionally, electrostatic interactions are modulated by the increased permittivity for silicon, which is ϵ≈ 11.68 ϵ_0. 
We utilize an analytic expression for the electrostatic potential in a quantum well from a square gate of radius a at a distance z between the gate and the well <cit.>





    V_i(x,y) = V_g [g(x-a, y-a) + g(x-a, a -y) 
    
     +  g(a-x, y-a) + g(a - x, a -y)]

where

    g(x,y) = 1/2πarctanxy/z√(x^2 + y^2 + z^2)

and V_g is the voltage applied to the gate electrode. An example of a possible gate architecture with its accompanying electrostatic potential is shown in Fig. <ref>.
For simplicity of demonstration, we will focus on square gates. 
For all of our calculations, the plunger gates have a radius of 25 nm and the barrier gate is 30×50 nm as shown in Fig. <ref>a.
We emphasize that these choices are typical for the currently experimentally used Si-based quantum dot qubits, and our qualitative results are independent of these specific numerical choices for the confining potential. 
Our approach can be implemented using any electrostatic potential.
Gates of different shapes, sizes, and pitches can be used, the only thing needed is a grid of the confining potential producing the dots.
Our choice of potential is arbitrary and could use any level of sophistication. 
The potential used in this technique could easily be solved by carrying out the standard Schrodinger-Poisson self-consistent calculation using the relevant lithographic geometry. 
However, we use the model potential confinement above for the dots instead of doing a numerical Schrodinger-Poisson approach because we are presenting a general theory instead of modeling specific devices.
We note as an aside that the Schrodinger-Poisson technique is essentially a one-particle problem, and is not particularly computationally demanding compared with the FCI technique we employ for extracting the Hubbard parameters as discussed below.

In the case of multi-electron dots, calculating the exact many-body wavefunction for the double quantum dot becomes prohibitively expensive (even for our model confinement potential). 
This is because the minimum number for each dot is three in order to have a valence electron, and N=6 between the two dots is expensive enough to make an FCI calculation impractical. 
This leads us instead to treat the intradot and interdot interactions separately. 
We assume that the correlation effects between the two dots are small enough to be safely ignored for the FCI calculations, and we calculate the FCI ground state for each dot separately. This is also consistent with the atomic orbital approach for calculating the Hubbard parameters since the generalized Hubbard approach to spin qubits implicitly assumes this individual dot independence we are using in our FCI theory.

However, the single dot potential in Eq. <ref> is incomplete for dot-specific FCI calculations because these single dot potentials are unaffected by changes to the barrier gate.
In reality, we expect a lowered barrier gate to significantly increase wavefunction overlap between the two dots and that a raised gate would decrease that same overlap.  
We ensure that each single dot potential is both localized and accurately reflects the full potential—from both plunger gates and the barrier gate—at that site by taking the full potential and rapidly tuning it to zero as it approaches the second dot. 

    V_eff(x,y) = V(x,y) exp[-max(0.0,(x - d/4))^2/(d/2)^2]

where 

    V(x,y) = V_L(x,y) - V_X(x,y) + V_R(x,y).

This approach is illustrated in Fig. <ref> and we refer to these adjusted dots as effective single dot potentials.





Configuration interaction (CI) methods require expressing the many-body Hamiltonian in the basis of Slater determinants.
Each Slater determinant represents a different possible electron configuration, and the Coulomb term in Eq. <ref> couples different configurations, hence the term configuration interaction.
Given M orbitals and N electrons, there are K = 2MN possible determinants.
The process is described as full configuration interaction when all K determinants are employed.
Given the fact that the number of determinants K grows very quickly, we carefully choose our orbitals to keep the computational cost within reason and avoid memory constraints.
Rather than use a predetermined basis of orbitals for all calculations, we use the eigenstates of the single-particle Hamiltonian in Eq. <ref>.
This approach has already been used effectively for spin qubit modeling of single electron dots by Anderson et al. <cit.>.
This allows us to minimize the number of orbitals M which must be kept.
Therefore, for each choice of gate voltages and pitch, we solve for the single-particle eigenstates of each dot using exact diagonalization in a truncated Fourier basis consisting of N_F basis functions.
The eigenstates are approximated well after about N_F = 500 basis functions, and the eigenenergy convergence is shown in Fig. <ref>a.



Including enough orbitals M is a critical part of a successful calculation, and we demonstrate the FCI ground state energy convergence in Fig. <ref>b.
Increasing M until the output is consistent and predictable is mandatory (and is, in fact, the definition of self-consistent FCI), as it is very difficult to know a priori how many orbitals will be necessary for sufficiently converged Hubbard parameters.
Ground state convergence is not by itself sufficient for converged Hubbard terms, however.
Separating out parts of a system's wavefunction and projecting to the Hubbard space requires higher-order terms arising from the excited states to be included. 
Early truncation can lead to erratic, numerically unstable values of the Hubbard parameters.
We have found keeping the first M=30 Slater determinants to be good practice for our choice of parameters, but this has to be ensured in each calculation since, depending on the details, more or fewer orbitals may be necessary for the computation.

Next, we need to extract the effective single-particle states from the FCI wavefunction.
For a single Slater determinant, this is simple.
However, such an approach would ignore correlational effects within the dot.
Separable approaches to the many-body system will never capture dynamical correlations, but important contributions from exchange and static correlation can still be incorporated.
Fortunately, when N=3, it becomes possible to isolate the valence electron well in such a way as to ensure that the two core electrons resemble a singlet state as much as possible.
This process consists of taking the one-particle reduced density matrix (1RDM) <cit.>, which is given as

    γ_ij = ⟨Ψ | â_i^†â_j | Ψ⟩,

where | Ψ⟩ is the ground state FCI wavefunction.

The eigenvectors of this density matrix represent the natural orbitals of the system, and the eigenvalues are their occupation numbers.
Because there are three electrons in the dot, one of them can be isolated by its spin, which without loss of generality, we classify as spin up.
The natural orbitals that are spin up, and their occupation numbers represent one of the core electrons, and we construct a density matrix to describe its state.
We then can use the remaining spin-down orbitals to construct another state that resembles the core electron as best as possible. 
The remaining state is that of the valence electron.
This is similar to the occupation number representation of a Slater determinant since it is not important which electron interacts with which.
Only a correct understanding of the states and their occupation is required.
Consequently, the many-body state is separated into three single-particle operators.
Each single-particle density operator's spin is also well-defined, which maps well to the Hubbard representation.









These single-particle density operators then need to be orthogonalized with respect to the same operators on a neighboring dot since they will have a nonzero overlap for any finite pitch.
This is usually done by performing a Löwdin orthogonalization process, which is symmetric with respect to its arguments.
Under the Löwdin procedure, the vectors are rotated so that no single vector is rotated more than the others.
This way, all new vectors will resemble their original starting points as much as possible. 
This is the standard orthogonalization routine in quantum chemistry <cit.>.
Löwdin orthogonalization differs from the more familiar Gram-Schmidt process, which is strongly asymmetric in its treatment of the initial vectors.
However, a direct application of the Löwdin procedure is obviously inappropriate, since density matrices do not form a vector space. 
Therefore, we again work with the natural orbitals of each density matrix.

To orthogonalize our set of single particle density matrices, we use the density operator Uhlmann fidelity F(ρ, σ) as the distance metric over which we define orthogonality.

    F(ρ,σ) = tr(√(√(σ)ρ√(σ)))

Fidelity between two density operators F(ρ,σ) = 0 if and only if ρ and σ have orthogonal support <cit.>.
Löwdin orthogonalization on the eigenvectors of the single-particle density operators produces density matrices of single-particle states on neighboring dots with zero fidelity.

This orthogonalization is done as follows: We perform a change of basis from the 2M basis orbitals (M from each dot), which are not mutually orthogonal, to orthogonalized versions of these basis states.
This is a passive transformation to place density matrices from neighboring dots on equal footing and to work in an orthonormal basis. 
We then find the natural orbitals that contribute most to each density matrix and orthogonalize them with respect to each other using the Löwdin procedure.
Orbitals are iteratively included until overlaps of states from neighboring dots reach a set threshold ε = 1×10^-10. 
This is done because as more (irrelevant) orbitals are included in the orthogonalization procedure, the resulting states resemble their original versions less and less.
This lets us directly extract the relevant Hubbard parameters for obtaining the charge stability diagram. 
We need only to modify the initial equations to handle density matrices rather than wavefunctions.

    F_ij = tr√(√(ρ_i)Ĥρ_j Ĥ√(ρ_i))


    G_ijkl = tr√(√(ρ_i ⊗ρ_j)Ĥρ_k ⊗ρ_l Ĥ√(ρ_i ⊗ρ_j))

where ρ_i corresponds to the state c^†_i |0⟩ in the Hubbard model.
With these quantities, we can efficiently simulate the charge stability diagram for a given gate architecture.





§ RESULTS




We take a one-dimensional slice in the charge stability diagram along the V_1 = V_2 diagonal and calculate the FCI ground state of each dot for N=1,2, and 3 to find the correct electron number. 
These results are shown in Fig. <ref>.
The charge transitions of V = 0.30 and V=0.58 mark the boundaries of the charge regimes which can be used for further calculations. These transition points will depend on gate layout and barrier gate strength but are relatively cheap to determine for N ≤ 3 electrons since the ground state energy converges quickly in M compared to Hubbard values.

We first calculate the tunnel coupling as a function of plunger gate voltage, barrier gate voltage, and dot pitch.
Fig. <ref> displays our results.
We find that tunnel coupling varies significantly over different orbital states and that t_3, which corresponds to accessing the first p-orbital of each dot, sees a  50x increase in tunnel coupling strength relative to the N=1 case for modest values of pitch and barrier gate.
This is sensible since the electrons are less localized to their original dot as they access higher orbital states, and we expect to see the overlap increase exponentially in such a regime. 
It is worth noticing that the value of t fluctuates significantly within each charge state, with a monotonically decreasing behavior as the plunger gates are increased.
This is physically logical, as we would expect the states to become more localized. 
However, for the purposes of modeling charge stability diagrams, and specifically the transitions between charge states, this effect is minor compared to the order-of-magnitude shifts that occur at the charge transitions.

The effects of the barrier gate are apparent, with increased values of V_X immediately diminishing the result of t_3 ≫ t_1. 
This behavior is only observable because of our effective single dot potentials, which take into account the barrier gate strength in calculating the center of the single dot well. 
As we increase the barrier voltage, it appears that tunnel couplings for the n=3 electron drop below that of the core electrons. 
This is likely only a numerical artifact, as the tunnel couplings in that regime are small enough that our separation and orthogonalization procedure may fail to capture the correct (but very small) value. 
Numerical errors in such exponentially small tunnel coupling values are not a problem for obtaining accurate charge stability diagrams.

Analyzing the tunnel couplings calculated as the dot pitch is varied (Fig. <ref>) also illustrates that this regime of very large t_valence is quickly suppressed for larger pitches or barriers. Here it appears that d = 90 nm is a key value (for our parameters), after which the tunnel coupling starts to drop precipitously. Similar calculations could be useful for estimating the spatial distribution of the valence electrons in individual devices.





The effect of these occupancy-dependent tunnel couplings is seen when generating the charge stability diagram.
In Fig. <ref>a, we calculate a charge stability diagram using only Hubbard parameters calculated from a double dot system with plunger gates at 0.2V. 
The Hubbard parameters are applied over the whole stability diagram.  
Switching on level-dependent tunneling leads to the picture in Fig. <ref>b, as the increased tunnel coupling leads to greater corner smoothing for higher charge states.
Because the major changes in t occur at charge transitions, we model t in the Hubbard model as being piecewise constant over each Coulomb diamond. This already represents a significant improvement over the widespread practice of assuming that the tunnel coupling is constant over the entire stability diagram. If a transition involves tunneling from one dot with n electrons into another with m electrons, the tunnel coupling corresponding to that interaction is t_j, where j = max(n,m+1). 
This occupancy-dependent tunneling is a new qualitative feature of our theory compared with the existing charge stability theories in the literature.

Because of the dramatically higher levels of tunnel coupling, points on the charge stability diagram that are one charge state in Fig. <ref>a could be another one entirely in Fig. <ref>b. 
Another feature in the level-dependent coupling diagram is that for N > 3, the individual Coulomb diamonds become less pronounced, and the diagonal slices corresponding to the double dot total charge become more straight and featureless (see Fig. <ref>). 
This feature is often seen in the experimental data <cit.>.
Experimental data usually only shows changes in the total charge of the system, though occasionally charge transitions at the dot level are also visible for a small dot occupation. Fig. <ref>b reflects how level-dependent tunneling may appear in total charge transition diagrams.




We performed a similar analysis for the onsite Coulomb term U in Fig. <ref> and the Coulomb coupling term V in Fig. <ref>.
The onsite Coulomb term U is unique in that it does not couple states of neighboring dots like t or V, and depends entirely on the properties of a single dot.
For that reason, it is calculated differently. Rather than using atomic orbitals to calculate the overlap of the charge distributions, we simply calculate the expectation value of the Coulomb Hamiltonian for our FCI single-dot ground state.
Fig. <ref> shows that the pitch and barrier gate have very small effects on U. However, as before with the tunnel coupling, we see significant fluctuations within a given charge state.
It is reasonable and even expected that the Coulomb energy increases as the dot becomes more localized.
However, in keeping with our purpose of modeling charge transitions, U is assumed to be constant over the charge stability diagram, since we do not see dramatic shifts in the average value of U between Coulomb diamonds, unlike the tunnel coupling.

Our calculations of the Coulomb coupling term V are shown in Fig. <ref>. 
The Coulomb energy between valence electrons of neighboring dots appears to have large shifts at charge transitions and to remain relatively constant within each regime.
As the wells become deeper within a charge regime, we do not expect V to change significantly, whereas we would expect larger changes to V when the dot wavefunctions have larger overlaps, such as when p-orbitals are occupied, or the barrier gate or pitch is lowered.
This is exactly what we see in the results, with p-orbitals having significantly larger Coulomb couplings.
We also observe V decreasing with barrier gate voltage and pitch in a manner generally similar to that of the tunnel coupling. 

It would appear that a piecewise constant approach to V would be appropriate, similar to our approach with the tunnel coupling. However, even if such an approach is pursued, the relevant quantity is the average value of V = 1/n_1 n_2∑_i,j V_ij. We employ this approach, and our charge stability diagrams use a constant value of V which is the average of the Coulomb couplings, since this average stays approximately the same across a single charge stability diagram.






The tunnel coupling varies dramatically throughout the stability diagram, and still only has a small (but nonetheless important) effect on the final stability diagram. Therefore, it is to be expected that the changes to the Coulomb terms, both the intradot U and the interdot V repulsion, can be safely ignored within a single stability diagram since the corrections in charge distribution that arise from accessing higher orbital states have only a minor effect on the Coulomb energies. Nonetheless, it is worth considering the effect that the barrier gate and dot pitch have on the entire charge stability diagram. We show these effects in Fig. <ref>.

The fact that V significantly increases as the barrier gate is lowered means that charge stability diagrams can look drastically different for the same gate geometry but with a slightly lower barrier gate voltage. Lowering barrier gate voltage also tends to decrease U slightly, but the increase in V is much more dramatic.  We again stress the importance of our effective single dot potentials for capturing these effects. Fig. <ref> clearly demonstrates that given a situation with enhanced t_3 coupling, this behavior can be suppressed by raising the barrier gate, and thereby returning the charge stability map to one where N=3 Coulomb diamonds more closely resemble those corresponding to N=1.


Our charge stability plots are generated in real space, starting from the (μ_1, μ_2) space and including neighboring dot occupation and electrostatic crosstalk through the interdot Coulomb effects with the standard formulas from the classical capacitance model.


    μ_1    = |e|[α_1 v_1 + (1-α_1) v_2] 
    μ_2    = |e| [(1-α_2) v_1 + α_2 v_2] 
    α_1    = (U_2 - V)U_1/(U_1 U_2 - V^2) 
    α_2    = (U_1 - V)U_2/(U_1 U_2 - V^2)


The alternative is to ignore these effects and plot the stability diagram using these theoretical quantities as the axes. 
These quantities are known as “virtualized”.
Gate virtualization results in “straightened out” stability plots with horizontal or vertical charge boundaries. 
It is more meaningful for theorists to express the stability diagram results in physical units to better compare their work to experimental results since the virtualization procedure distorts the stability diagram in such a way that the effects of modeling choices can be exaggerated. 
It could also be beneficial for experimentalists to provide virtualized charge stability diagrams to better connect with theory. An example of a virtualized stability diagram is shown in Fig. <ref>.










§ CONCLUSIONS

In conclusion, we have presented a multi-band and multi-electron interacting model of gate-defined double quantum dots and simulated charge stability diagrams for these systems using realistic electrostatic potentials, full configuration interaction, and effective Hubbard model mapping.
Our approach is compatible with confining potentials at any level of sophistication, including potentials obtained from Poisson-Schrodinger solvers as the starting point (replacing our model analytical confining potentials).
We separated the double dot system using an effective single dot approach that accounts for the barrier gate voltage.
We then utilized full configuration interaction to account for many-body interactions within each quantum dot.
Single particle density matrices were computed using natural orbitals of the 1RDM, and Hubbard couplings were calculated as the Uhlmann fidelities of the corresponding operators.
These parameters of the generalized Hubbard model enabled us to predict large regions of the charge stability diagram with improved accuracy compared to classical capacitance models or single-band Hubbard models. 

Our results highlighted the significant enhancement of tunnel couplings for valence electrons in dots with three electrons compared to single-electron dots. 
This difference was evident in the charge stability diagram for higher occupancy Coulomb diamonds. 
We extrapolate this trend to a large charge stability diagram in Fig. <ref>.
As dot occupancies increase, so do the tunnel couplings, and total charge regions become increasingly featureless and eventually become almost parallel lines.
Unlike the tunnel coupling, the average values of the onsite Coulomb term U and the interdot Coulomb coupling V were not found to vary significantly within the charge stability diagram.
We also found that increasing barrier strength or dot pitch reduces the effect of increased tunnel coupling, and this trend allows researchers to make their Coulomb diamonds as periodic as possible so that they see similar behavior using N=3 dots as they do N=1 dots.

Our FCI approach is only applicable in regimes where it is appropriate to treat the dots as separate subsystems, where wavefunction overlap is small. 
Double dot systems with large detunings or negative barrier gates are not compatible with this technique.
Likewise, calculations of interdot exchange interaction rely on extremely precise pictures of the correlational energy of the entire double dot system, and cannot be performed in this formalism, even when tunnel couplings are small. 
Although calculating the FCI energy of the double dot system as a whole is an excellent approach for such a problem <cit.>, this would likely be computationally feasible only for N<3 electrons in each dot.

This work opens the door to many exciting future research directions, such as how background disorder in the form of charge noise or magnetic impurities impacts the Hubbard parameters.
Evaluations of the resulting charge stability diagrams could be done by calculating the ground state charge states of the double dot system explicitly using FCI or another quantum chemistry approach. Although FCI would yield an extremely accurate answer, it probably could only be used for the first few Coulomb diamonds because of computational cost.
One advantage of this approach is that because it separates the dots, it could feasibly be applied to spin chains longer than two to investigate how next-nearest and next-next-nearest neighbor tunnel couplings scale compared to nearest neighbor couplings.

Our work demonstrates the importance of considering multi-electron effects in double quantum dot systems and the potential for the generalized Hubbard model to capture these effects. 
By accurately predicting the charge boundaries of the stability diagram, our approach provides valuable insights for tuning multidot devices and achieving specific electron number states. 
We believe that this technique will allow for cost-effective simulation of very high-quality charge stability diagrams, which should be relevant for constructing gating sequences, prototyping device designs, and generating data for machine learning algorithms. 
Overall, our findings significantly contribute to the understanding of charge stability diagrams of quantum dot semiconductor spin qubits and enhance their potential for scalable quantum computing applications. 



§ ACKNOWLEDGEMENT

This work was supported by the Laboratory for Physical Sciences.


