Interacting electrons on a quantum ring: exact and variational approach

We study a system of interacting electrons on a one-dimensional (1D) quantum ring using exact diagonalization (ED) and the variational quantum Monte Carlo (VMC) method. We examine the accuracy of the Slater–Jastrow-type many-body wavefunction and compare energies and pair distribution functions obtained from the two approaches. Our results show that this wavefunction captures most correlation effects. We then study the smooth transition to a regime where the electrons localize in the rotating frame, which for the ultra-thin quantum ring system happens at quite high electron density.


3
Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT Fogler and Pivovarov [14] and show that they are in good agreement. We discuss the localization of electrons in relative coordinates and find a smooth transition starting already at 1 < r s < 2.
This paper is divided as follows. In section 2, we discuss our model of the quantum ring system and go in considerable detail through the two methods used to calculate its ground state properties, ED and VMC. Then the results are presented in section 3. Our conclusions follow in section 4.

Model
We model the semiconductor quantum rings in the effective mass approximation by a many-body Hamiltonian where θ i is the angular coordinate of the ith electron and V(θ ij ) is the interaction potential between each pair of electrons. We assume that the ring is extremely narrow so that the electrons can only move at a certain radius, leading to a purely 1D ring. The Coulomb interaction for a ring system is written as where R is the radius of the ring, V C controls the strength of the interaction, and µ is a small parameter that eliminates the singularity at θ ij = 0. As the interaction between particles at the same position on the ring is V C /µR, one can think the regularization of the potential to give the ring a width of µR. The single-particle eigenstates of the noninteracting case are plane waves with energies where l = 0, ±1, ±2, . . . is the angular momentum quantum number. We scale the units so that the energy is measured in units of E 0 = h 2 /(2m * π 2 R 2 ) and length in units of R. The only parameter which is tuned is the interaction strength, and for the Coulomb interaction this can be translated into a change of the ring radius by the formula where a * B is the effective Bohr radius. The results presented in this paper are for interaction strength in the range of 1 to 10, which corresponds to a radius of 21 to 209 nm assuming GaAs 4 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT material parameters. For comparison, the quantum ring system studied by Keyser et al [1,15] had a radius of around 150 nm, which corresponds to an interaction strength of V C = 7.2 in our units. The width of their ring can be estimated to be around 20 nm, which is in reasonable agreement with the value of µ = 0.1 used in this study.
In the following, we will occasionally use the 1D r s -parameter, which is the average distance between electrons. In our model, it is directly proportional to the radius of the ring and thus also to V C r s a *

ED
The ED method, originally called configuration interaction, is a systematic scheme to expand the many-particle wavefunction using the noninteracting states as a basis. This method traces back to the early days of quantum mechanics, to the work of Hylleraas [16]. The calculation of matrix elements for the many-body case was originally derived by Slater [17,18] and Condon [19] and developed further by Löwdin [20]. The use of the term ED for this can be seen as an attempt to replace the term configuration interaction with a normal mathematical term [21]. This term might in some cases be misleading, as truly exact results are obtained only in the limit of an infinite basis. As a simple example to start from, consider a single-particle Hamiltonian split into two parts H = H 0 + H I , where the Schrödinger equation of the first part is solvable and the wavefunctions φ i (r) form a orthonormal basis. The solution of the full Schrödinger equation can be expanded in this basis as ψ(r) = i α i φ i (r). Inserting this into the Schrödinger equation results in Using (7), multiplying with φ * j and integrating gives for every j. This can be written as a matrix equation where H 0 is a diagonal matrix whose ith element is ε i , and the jith element of H I is φ * j (r)H I φ i (r) dr. The vector α contains the values α i . In this way, the Schrödinger equation In a similar fashion, one can split the many-particle Hamiltonian into two parts, typically a noninteracting and an interacting one. The solution to the noninteracting many-fermion problem is a Slater determinant formed from the N lowest energy eigenstates of the single-particle Hamiltonian. The basis typically chosen for the interacting problem is a Slater determinant basis constructed from M single-particle states. For an N-particle problem, M N different manyparticle configurations can be constructed by occupying any of the M available states. Each of these determinants is labelled by the N indices of the occupied single-particle states. Due to interactions, other configurations than the one of the noninteracting ground state have a finite weight in the expansion of the many-particle wavefunction.
The Schrödinger equation maps to matrix form in a similar way as in the preliminary example above. The calculation of the matrix elements of the Hamiltonian are simpler using second quantization and the occupation number representation. We will only state the results here.
The noninteracting Hamiltonian matrix H 0 in the determinant basis is diagonal with the ith element equal to E i = N j=1 ε i j , where i j labels the occupied single-particle state. The elements of the interaction Hamiltonian matrix can be found by a straightforward calculation using the anti-commutation rules of the creation and annihilation operators. The only nonzero matrix elements are between configurations that differ at most by two single-particle occupations. If both of the two configurations have two single-particle states occupied that are unoccupied in the other one, the interaction matrix element is where p and q are single-particle states occupied in the left configuration and unoccupied in the right one, and similarly for s and r with right and left interchanged. The sign of the matrix element depends on the total number of occupations between p and r and between q and s due to anti-commutation of the creation and annihilation operators. If this number is odd, the sign is minus, otherwise plus. For the Coulomb interaction and the delta functions result from the spin part summation. In the case of a 1D quantum ring this matrix element has a closed form where Q(x) is the Legendre function which can be written in terms of Gauss's hypergeometric function.
If the two configurations have a difference of only one occupation, the final result for the interaction matrix elements is where q is occupied in the left configuration and unoccupied in the right one, and for r the other way around. The sum over i is over orbitals that are occupied in both configurations. Again the sign of the matrix element depends on number of occupations between the differing orbitals.
It is easy to see that the matrix elements in (15) vanish for a rotationally symmetric system, as a change in only one occupation cannot conserve the total angular momentum. The diagonal elements of the interaction Hamiltonian matrix are simply where the sum runs over the occupied states.
To summarize, the noninteracting part results in a Hamiltonian matrix which is diagonal, and the interaction Hamiltonian couples configurations that differ by at most two occupations. Now that we know the Hamiltonian matrix elements, we can sketch the basic ED procedure as follows: First solve for M eigenstates of the single-particle Hamiltonian H 0 . After those have been found, calculate the two-particle matrix elements V ijkl of (13). Construct the M N N-particle configurations from M single-particle states. Here, configurations with wrong symmetry can be rejected. For example, one can pick only configurations with certain z-component of the total spin or some other good quantum number, like the angular momentum in our case. To construct the Hamiltonian matrix, calculate the interaction matrix elements of the Hamiltonian H I between the configurations. Construct the diagonal noninteracting Hamiltonian matrix elements from the single-particle eigenvalues. Finally, diagonalize the Hamiltonian matrix. The main problem of the ED method is the exponential computational scaling of the method as a function of the number of particles. An other problem is the convergence rate as a function of basis size.
To compare the ED and VMC methods for the ring, we study the total energy and the pair distribution function of the electrons in the ring. The pair distribution function of the system is defined by where n σ (r) is the spin σ particle density and † σ (r) and σ (r) are the field operators. It describes the probability that, given an electron of spin σ at r, an electron of spin σ is found at r . It therefore gives an idea of the strength of correlation between electrons. The pair distribution function can be derived in a similar way as the Coulomb matrix elements using the anti-commutation rules of the creation and annihilation operators.

VMC method
The QMC methods are among the most accurate ones for tackling a problem of interacting quantum particles [22]. Often the simplest QMC method, namely, the VMC, is able to reveal the most important correlation effects. In many quantum systems further accuracy in, e.g., the energy is needed, and in these cases methods such as the diffusion QMC allow one to obtain more accurate estimates for various observables.
The VMC strategy of solving the ground state of an interacting system of N particles (defined by the Hamiltonian operator H and the particle statistics) starts by constructing the variational many-body wavefunction . Then, all the observables of the system can be computed from it. For example, the energy E, which is higher than the ground state energy E 0 , can be obtained where R is a vector containing all the coordinates of the N particles. If the particles have d degrees of freedom, the integrals here are N × d dimensional. As one would like to have N reasonably large, the integrals are typically so high-dimensional that these are most efficiently calculated using a Monte Carlo strategy. This is based on the limit where the points R i are uniformly distributed throughout the volume V . The error in the evaluation of the integral decays as ∝σ f / √ N R , where σ f is the standard deviation of the function f . A typical trick to reduce the error is to make σ f smaller by dividing the function f by a similar function g (that we assume to integrate to one). Then, one can rewrite the sum as where now the points R g i are distributed as g. The error in this approach is now smaller, as the function g was chosen such that the standard deviation of f/g is smaller than that of f . Now in VMC, one typically generates the random set of N R coordinates {R i } N R i=1 so that they are distributed according to | | 2 , and one obtains for the energy where we have defined the local energy E L to be One should note that E L is a function of the coordinates, and it defines the energy at this point in space. If the wavefunction solves the Schrödinger equation, E L equals the true energy of the system. From the Monte Carlo point of view, the error in determining the energy is related to the standard deviation of the local energy, and for this reason the statistical error in the VMC energy decreases as the wavefunction becomes more accurate. The sampling of the random points R i can be done by, e.g., using the basic Metropolis algorithm. In that, the ratio of the sampling function at the ith and the previous step: | (R i )| 2 /| (R i−1 )| 2 , needs to be calculated and the trial step is accepted if the ratio is larger than a uniformly distributed random number between zero and one.
A typical, and also the most commonly used VMC trial many-body wavefunction, is the one of a Slater-Jastrow type: where the two first factors are Slater determinants for the two spin types, and where j(θ ij ) is a Jastrow two-body correlation factor. The determinants contain N single-particle orbitals, where the ones in different spin determinants can be the same. These determinants solve in some cases a noninteracting or a mean-field Schrödinger equation, but in general, one can choose these orbitals rather freely. This form of a wavefunction has proved to be very accurate in many cases [22]. In quantum dots, the Jastrow-Slater wavefunction captures most correlation effects very accurately [9]. The only possible exceptions are found in the fractional quantum Hall effect regime, where only some of the states can be well approximated by the Jastrow correlation.
Interesting examples of successful cases are the Laughlin states [23]. One can generalize the Jastrow-Slater wavefunction to contain more than one determinant per spin type. This leads to an increased computational complexity, and for this reason it is often avoided. Another generalization can be to replace the single-particle coordinates in the determinants by a set of collective coordinates, e.g., . This also leads to complications, this time in the calculation of the ratio of the wavefunctions in the Metropolis algorithm.
The determinants in our VMC trial wavefunction are constructed from the noninteracting single-particle states, and for the two-body Jastrow factor we use here a simple form of where α are variational parameters, different for electron pairs of same and opposite spin type. An important ingredient in VMC is the optimization of the variational parameters of the wavefunction [9]. For an efficient optimization of the variational parameters, we need the derivative of the energy with respect to these parameters. Our wavefunction is complex-valued, so we need to generalize the result for the real case [24] to read: the probability distribution | | 2 . This notation is also used in the following. In our wavefunction, the parameters α i are inside an exponential, real-valued Jastrow factor: where ≡ (θ 1 , . . . , θ N ), α ≡ (α 1 , . . . , α M ). As a result we have Using (28), we can rewrite (26) as

Convergence and configurations
Before presenting the actual results, we first present an analysis of the accuracy of the ED method used. The only approximation done in our ED calculations is that we have a finite number of determinants in our expansion of the many-electron wavefunction. The actual number used is finally limited by the available computer resources. To give an idea of the accuracy of ED, we show in figure 1  which is about 0.00034% of the total energy. We can therefore be confident that our ED calculations have sufficient accuracy to be used to test the quality of the VMC approach.
When the electrons interact weakly, the many body state has one major configuration, that of the equivalent noninteracting system. In this limit, the single-determinant Hartree-Fock (HF) method gives reasonable results. The electrons arrange so as to minimize the total kinetic energy and the single-particle states are occupied as compactly as possible around the l = 0 state. The resulting configurations are shown in the upper part of figure 2. Each box in this figure represents one single-particle state, having angular momentum from l = 0 at the bottom, up through l = ±1, ±2, . . . , with increasing energy. For four electrons, the second shell formed by the angular momentum l = ±1 states is half-filled, and the exchange energy favours spinpolarization of the two electrons in that shell. When the interaction strength grows, the relative importance of the kinetic energy decreases with respect to the interaction energy. The most important configuration of the many-body wavefunction of strongly interacting electrons has closed shells, and thus zero angular momentum, for each spin type. The configurations are then either a half-filled shell of one spin type or a closed shell. The resulting configurations are shown in the lower half of figure 2. For N = 4 and 6 electrons, the state is already a closed-shell at weak interaction and it remains the major configuration in the limit of strong interaction. The odd N cases have open-shells at weak interaction and make a transition to closed-shell configuration at some point.

VMC accuracy
The total energy of six electrons on the ring is shown in figure 3. We have done the calculations for the L = 0 state, and a total spin S = 0 (groundstate) and S = 2, and show the results obtained from ED, VMC and HF as a function of the strength of interaction (proportional to the ring radius). Over the whole range, ED and VMC results agree very well, although the agreement improves slightly with increasing V C . The HF results are reasonably accurate only at the lowest interaction strength and already at V C = 2 (corresponding to R = 4 a * B and r s = 2.1 a * B ) shows considerable deviation from the ED results.
To compare in more detail the results of ED and VMC, a convenient measure of the quality of the trial wavefunction is the percentage of correlation energy captured by VMC, defined as where E HF is the HF energy. In figure 4, we show the results for the L = 0 state of six electrons. At all interaction strengths, VMC captures most of the correlation energy. It is more accurate for the S = 2 state than the S = 0 state, although the difference decreases with increasing strength of interaction. One would expect VMC to become less accurate for strong interactions, because the construction of the VMC wavefunction starts from a single configuration which is later multiplied by a Jastrow factor, whereas the many-body wavefunction contains multiple determinants and their weights get more evenly distributed with increasing V C . This is not the case, and at the highest V C the deviation in the VMC groundstate energy lacks such a single-configuration limitation. For other particle numbers the results are shown in table 1 for low and high interaction strength, V C = 1 and 10. The VMC results are very accurate for the particle numbers we have considered. In the worst case, VMC captures 96% of the correlation energy and the total energy is roughly 0.7% higher than that found by ED. The accuracy improves with increasing interaction strength and is    (7) better for lower particle numbers. This trend is broken by cases where the major configuration has an open shell, which is always the case for odd particle numbers at weak interaction. This is demonstrated by five electrons, where VMC captures slightly less of the correlation energy than for six electrons. In figure 5, the total pair distribution function of an N = 6 electron ring is shown at weak and strong interaction. As in the case of the groundstate energy, the agreement between VMC and ED is good. The left panel shows clearly how the total pair distribution functions of different spin states become identical when the electrons interact strongly. This is one indication of spin-charge separation, when the spin and charge of a 1D system decouple.
The spin channels in the pair distribution function are shown in figure 6. For the S = 2 state, VMC results are almost identical to those from ED. In the case of S = 0, however, the internal structure is inaccurate at weak interaction and completely wrong at strong interaction. This is further shown in figure 7 where we have plotted which for an antiferromagnetic alignment of spins should give equidistant peaks (up spins) and troughs (down spins). This is indeed the case, and we furthermore see a decay in spin correlations with increasing interelectron distance. The VMC results show no antiferromagnetic structure. We believe the reason for VMC not being able to handle antiferromagnetic spin coupling is due to three-particle correlation which is missing in the VMC wavefunction. Despite VMC's shortcomings in accurately describing the spin correlations for an antiferromagnetic spin structure, the total pair distribution function and the groundstate energy are correctly described.

Spin and charge in the Wigner crystal regime
Like the pair distribution function, the total energy of the S = 0 and 2 states becomes identical in the limit of strong interaction. In figure 8, we show the energy difference of these states for six electrons, as a function of V C found by ED, VMC and HF. The inset shows the energy difference found by ED on a logarithmic scale. As expected, HF fails seriously already at low values of V C and predicts a groundstate of nonzero spin-polarization at V C just below 1 (r s ≈ 1). VMC follows the ED results more closely, but at V C = 4 predicts a change in spin. ED results indicate, however, that the groundstate spin remains S = 0 over the whole range of V C although the energy difference approaches zero. Interestingly, the difference between the ED and VMC results remains approximately constant for V C = 4 . . . 10. Figure 8 shows that in the strong interaction regime, the energy scale of spin dynamics becomes very small when compared to the energy scale of the orbital motion. This is one indication of a strong spin-charge separation in the system. It is further manifested in the pair distribution function, as shown in figure 5, which is the same for both spin states, S = 0 and 2, at stronger interaction. The spin part can be described by the antiferromagnetic Heisenberg model [13] with spin interaction with a nearest-neighbour coupling constant J. In a recent paper, Fogler and Pivovarov [14] derived a relation between the coupling constant of thin quantum rings and the pair distribution function of a spin-polarized 1D system g 0 (x) where R * is the thickness of the ring, in our case µR. For ultrathin wires exhibiting strong spin-charge separation, this relation holds even at moderate electron densities, r s ∼ 1. They furthermore derive an expression for the pair distribution function, which leads to where κ and η are constants. The coupling constant J is proportional to the difference in energy of the S = 0 and 2 states, E. In our model, r s = 2πV C /N, so we can write In the inset of figure 8, we have plotted the logarithm of E as a function of √ V C . Our results agree very well with the results of Fogler and Pivovarov [14], as the curve becomes linear even at low values of V C , V C ≈ 1.7 (r s ≈ 1.8). From our data we can estimate the slope, and find that η ≈ 2.5 which is reasonably close to the value 2.7979 found by Fogler and Pivovarov [14]. However, it must be noted that our results are for a fixed value of µ = 0.1, so that the thickness of the ring scales linearly with increasing interaction strength as The pair distribution function of six electrons of zero total angular momentum and spin is shown in figure 9 for several values of V C . At V C = 1, the distribution is rather flat, but already at V C = 2 there is a clear peak structure indicating partial localization in relative coordinates. As the system is finite, a phase transition to a Wigner crystal does not take place. However, there is a smooth transition to a regime where a Wigner-crystal-like behaviour is found. The system is translationally invariant and the electron density remains uniform, the localization occurs in the relative frame. In 2D quantum dots, the electrons are found to be localized at [25] r s = 4, and the transition is preceded by a change in the groundstate spin [26]. A QMC study of the 3D electron gas showed an abrupt change in the evolution of the correlation energy with density [27]. We do not see such a change in our system. One should also note that the 3D QMC work has different wavefunctions for the liquid and solid phase, whereas we use the same wavefunction in both cases.
Further insight is obtained from the value of the pair distribution function at θ = 0, which is the contact probability. At low values of the strength of interaction V C , it has a nonzero value, as electrons of opposite spin are not restricted by the Pauli exclusion principle. As V C grows, the strong Coulomb repulsion reduces g ↑↓ (0) to the point where the electrons become localized between their neighbours on either side. A closer look at the contact probability reveals an indication that the quantum ring system undergoes a smooth transition at densities between 1 < r s < 2. In figure 10, the logarithm of the contact probability is shown as a function of the interaction strength on a square root scale. For V C 2, the linear fit indicates that g ↑↓ (0) falls off as exp(−ζ √ V C ), where ζ is a constant. At low interaction strength, the behaviour is closer to exp(−ζ V C ), the parabolic curve shown. A related behaviour has been found for a different quantity in quantum dots [25].

Conclusions
We have studied interacting electrons in a 1D quantum ring. We compared the results obtained from VMC with those of ED. The variational wavefunction was taken to be of a simple Slater-Jastrow form: a single determinant, the ground state of the noninteracting system, multiplied by a two-body correlation factor. The results obtained indicate that this simple wavefunction is able to capture the important correlation effects, and in most cases even with a high accuracy. The major source of error in VMC is the spin energy scale and spin correlation. The spin energy scale, however, becomes very small already at rather weak interaction and the error in the groundstate energy is correspondingly small and decreases with increasing V C . For an antiferromagnetic spin correlation, VMC does not capture the internal spin structure, but the total pair distribution function is accurate. From this, one can conclude that the Jastrow factor is very efficient in capturing the total energy and essential correlation effects induced by the interaction between electrons. One might question the importance of the determinant part of the VMC wavefunction. We believe that having a good choice for this is crucial for the success of the VMC simulation. Even if the weight of this single configuration is not large in the exact expansion when the interaction is strong, it specifies the symmetry of the VMC many-body state. One should note that this is also the case for quantum dots, where the noninteracting singleparticle states can be used in the determinant part of the VMC wavefunction, and they are even the optimal ones in the case of a closed-shell configuration [9]. The reason, both in quantum rings and dots, is the high symmetry of the problem: in a parabolic dot, the centre-of-mass and relative motion decouple, also in the many-body wavefunction. Thus if one starts from a noninteracting state that is a ground state of the centre-of-mass motion, the interactions only change the relative motion part of the many-body wavefunction. This means that if a Slater-Jastrow wavefunction is used, the Slater determinant part is not changed by the interaction between electrons. The case of a quantum ring is similar in the sense that the centre-of-mass and relative motion also decouple, but now because of the translational invariance of the system.
For a moderate strength of interaction, the spin correlation is well described by an antiferromagnetic Heisenberg model with coupling constant J. We showed that our results agree well with the expression for J derived by Fogler and Pivovarov [14] at r s 2. Their result is applicable already at r s ∼ 1 if a strong spin-charge separation is present in the system. This indicates that for an ultrathin quantum ring this separation does take place at a quite high electron density. The contact probability also shows that there is a change in the properties of the system starting at 1 < r s < 2.
In conclusion, a simple variational wavefunction, combined with Monte Carlo techniques to calculate the optimal wavefunction and observables of the system, results in an efficient computational strategy for interacting electrons in purely 1D quantum rings. Since a similar conclusion can be drawn for 2D quantum dots [9], we expect this to hold also for quantum rings of a finite width.