Orbital ferromagnetism in interacting few-electron dots with strong spin-orbit coupling

We study the ground state of $N$ weakly interacting electrons (with $N\le 10$) in a two-dimensional parabolic quantum dot with strong Rashba spin-orbit coupling. Using dimensionless parameters for the Coulomb interaction, $\lambda\lesssim 1$, and the Rashba coupling, $\alpha\gg 1$, the low-energy physics is characterized by an almost flat single-particle dispersion. From an analytical approach for $\alpha\to \infty$ and $N=2$, and from numerical exact diagonalization and Hartree-Fock calculations, we find a transition from a conventional unmagnetized ground state (for $\lambda<\lambda_c$) to an orbital ferromagnet (for $\lambda>\lambda_c$), with a large magnetization and a circulating charge current. We show that the critical interaction strength, $\lambda_c=\lambda_c(\alpha,N)$, vanishes in the limit $\alpha\to \infty$.


I. INTRODUCTION
The electronic properties of few-electron quantum dots in semiconductor nanostructures have been widely studied over the past decades [1][2][3]. Typically, electrons in the two-dimensional (2D) electron gas formed at the interface between different semiconductor layers are confined to a localized region in space by means of electrostatic trapping. The resulting confinement is usually well approximated by a parabolic potential with oscillator frequency ω, suggesting a simple 2D oscillator spectrum. However, Coulomb interactions are important in such devices, and their impact can readily be seen in transport spectroscopy [2]. Apart from the ubiquitous Coulomb charging effects, they are also predicted to induce a transition to a finite-size Wigner crystal of N electrons, the "Wigner molecule" [4,5], where the electrostatic repulsion suppresses quantum fluctuations and inter-electron distances are maximized [6][7][8]. The ratio between the confinement scale, l T = /m e ω, with the effective mass m e , and the Bohr radius, a B = 2 ε 0 /m e e 2 , defines a dimensionless interaction strength parameter [3], (1.1) Interactions are here described by the standard Coulomb potential, V (r) = e 2 /ε 0 r, where the dielectric constant ε 0 accounts for static external screening. The crossover from the weakly interacting Fermi liquid phase (realized for λ 1) to the Wigner molecule then happens around λ ≈ 1 and is known to be rather sharp despite of the finite-size geometry [5]. Due to the confinementinduced reduction of quantum fluctuations, the corresponding electron densities near the transition are much higher than the one required for bulk Wigner crystal formation [3].
Another modification of the 2D oscillator spectrum is caused by spin-orbit coupling. We here focus on the Rashba term caused by interface electric fields, which often is the dominant spin-orbit coupling and can be tuned by gate voltages [9]. Other types of spin-orbit coupling are expected to generate similar physics as described below, assuming that one can reach the corresponding strong-coupling regime. In particular, the model studied below applies directly to the case of Dresselhaus spinorbit coupling [9]. With the Rashba wavenumber k 0 , it is convenient to employ a dimensionless Rashba coupling, (1. 2) The single-particle spectrum of a dot with weak Rashba coupling, α 1, has been discussed, e.g., in Refs. [10,11].
In this paper, we study interacting few-electron quantum dots in the regime of large Rashba spin-orbit coupling, α 1. This regime appears to be within close experimental reach [22][23][24][25][26][27][28][29], and is also of considerable fundamental interest. In fact, many materials with strong spin-orbit coupling are known to realize a topological insulator phase [30,31]. Near the boundary of a noninteracting 2D topological insulator with time reversal symmetry (TRS), an odd number of gapless one-dimensional (1D) helical edge states must be present [30,31], where the spin is tied to the momentum of the electron. As we do not address magnetic field effects here, the Hamiltonian below enjoys TRS. Moreover, it is characterized by strong spin-orbit coupling, and it resembles a topological insulator in the absence of interactions.
Given the above developments, it is not surprising that several theoretical works [32][33][34][35][36] have already addressed the physics of noninteracting electrons in quantum dots with α 1. In this limit, the low-energy spectrum of a parabolic dot is well described by a sequence of almost flat Landau-like bands (see Sec. II A), E J,n ω n + 1 2 + J 2 2α 2 , (1. 3) with half-integer total angular momentum J and the band index n = 0, 1, 2, . . ., such that states with the same n but different J are almost degenerate. Equation (1.3) reflects the spectrum of a 1D (radial) oscillator plus a decoupled rotor with large moment of inertia. Assuming that the Fermi energy is within the n = 0 band, with corresponding Fermi angular momentum J F , the Kramers pair with J = ±J F has eigenfunctions localized near the "edge" of the dot. In fact, those states have the largest distance from the dot center among all occupied states, and form a helical edge with opposite spin orientation of the counterpropagating ±J F states [32]. By virtue of the bulk-boundary correspondence [30], the authors of Ref. [32] argued that a noninteracting dot with α 1 has features similar to the finite-size version of a 2D topological insulator. Indeed, time reversal invariant singleparticle perturbations, e.g., representing the effects of elastic disorder, are predicted not to mix opposite-spin states, and the helical edge is therefore protected against such sources of backscattering. In the finite-size dot geometry, however, the Z 2 invariant commonly employed to classify the topological insulator phase is not well defined.
For a dot with α 1, since the noninteracting spectrum is almost flat, one can expect that interactions have a profound effect. For instance, in lattice models hosting a topological insulator phase for weak interactions, Mott insulator or spin liquid phases emerge for strong interactions [37]; for the case of interacting bosons, see Refs. [38][39][40][41][42]. Moreover, the conspiracy of a single-particle potential with sufficiently strong Coulomb interactions can induce two-particle Umklapp processes destroying the helical edge state [43,44]. Motivated by these developments, we here study the ground state of interacting electrons in a quantum dot with strong Rashba spin-orbit coupling. We find it quite remarkable that the relatively simple Hamiltonian below captures such diverse behaviors as Wigner molecule formation, the presence of helical edge states, and -as we shall argue -the molecular equivalent of an orbital ferromagnet. This Hamiltonian is also expected to accurately describe semiconductor experiments, where recent progress holds promise of reaching the ultra-strong Rashba coupling regime. Let us now briefly summarize our main results, along with a description of the structure of the paper.
In Sec. II A, we present the single-particle model for the quantum dot, and summarize its solution for large Rashba parameter α. While our general conclusions hold for arbitrary radially symmetric confinement, quantitative results are provided for the most important case of a parabolic trap. We introduce a single-band approximation valid for weak-to-intermediate interaction strength, λ 1, and energy scales below ω, which allows one to make significant analytical progress. In Sec. II B, we then discuss the general properties of Coulomb matrix elements. The limit of ultra-strong Rashba coupling, α → ∞, is addressed in Sec. II C, where a simple analytical result for the Coulomb matrix elements is derived. For the resulting α → ∞ model, H ∞ , already weak interactions induce strongly correlated phases. The Coulomb matrix elements not included in H ∞ , arising for large but finite α, are addressed in detail in Sec. II D.
Next, in Sec. III, we present the exact ground-state solution of H ∞ for two electrons (N = 2). While the above discussion may suggest that a Wigner molecule will be formed, we find an orbital ferromagnetic state. The N = 2 ground state of H ∞ , see Sec. III A, is shown to be highly degenerate in Sec. III B. However, perturbative inclusion of Coulomb corrections beyond H ∞ , see Sec. III C, breaks the degeneracy and suggest the possibility of spontaneously broken TRS in an interacting N = 2 dot (for a more precise characterization of this phenomenon, see Sec. III), with a large value of the total angular momentum found already for weak interactions. The emergence of a finite magnetization [45], M s = 0, suggests a finite-size ("molecular") version of an orbital ferromagnet. This remarkable behavior appears at arbitrarily weak (but finite) interaction strength, with giant values of the magnetization. We estimate M s ≈ (λα) 1/4 , see Sec. III C. This highlights that the orbital angular momentum is behind this phenomenon, see also Ref. [46].
In Sec. IV A, we then present exact diagonalization results for the ground-state energy of N = 2 and N = 3 electrons in the dot for α = 10 and α = 15, going beyond the α → ∞ model H ∞ . We now find that only above a critical interaction strength, λ > λ c (α, N ), the dot develops a magnetization, M s = 0. The parameter λ c becomes smaller with increasing α, which is consistent with λ c (α → ∞) → 0 as obtained from H ∞ in Sec. III. In Sec. IV B, we then discuss Hartree-Fock (HF) results for particle numbers up to N = 10, where exact diagonalization becomes computationally too expensive. The HF results show qualitatively the same effects, indicating that orbital ferromagnetism represents the generic behavior of weakly interacting electrons in quantum dots with ultra-strong Rashba coupling. Finally, we conclude in Sec. V, where we also discuss perspectives for experiments. Additional details about the α → ∞ limit are given in an Appendix.

II. COULOMB INTERACTIONS IN A RASHBA DOT
A. Single particle problem We consider electrons in a 2D quantum dot with parabolic confinement in the xy plane. Including the Rashba spin-orbit coupling, the single-particle Hamilto-nian reads [9] where k = −i(∂ x , ∂ y ), r = (x, y), ω is the trap frequency (defined in the absence of spin-orbit coupling), and the positive wavenumber k 0 determines the Rashba coupling. With Pauli matrices σ x,y,z referring to the electronic spin, the Hermitian helicity operator, P h = (k y σ x − k x σ y )/k, has the eigenvalues ±1. In the absence of the trap (ω = 0), helicity and momentum are conserved quantities. Writing k = k(cos φ, sin φ), it is a simple exercise to obtain the P h -eigenspinors, Φ ± (φ), with conserved helicity ±1. The dispersion relation is then (up to a constant shift) given by 2 (k ∓ k 0 ) 2 /2m e . Low-energy states have positive helicity with and for given k ≈ k 0 , a U(1) degeneracy is realized, corresponding to a ring in momentum space.
In the presence of the trap, however, helicity and momentum are not conserved anymore. The system now has two characteristic length scales, namely the confinement scale, l T = /m e ω, and the spin-orbit length, 1/k 0 . Their ratio determines the dimensionless Rashba parameter α in Eq. (1.2). In this paper, we discuss the case α 1, where positive helicity states are separated from P h = −1 states by a huge gap of order 2 k 2 0 /m e = α 2 ω. As a consequence, negative helicity states can safely be projected away. Noting that the total angular momentum operator, J z = −i ∂ φ + σ z /2, is conserved, with eigenvalues J (half-integer J), the low-energy eigenstates of H dot for α 1 have the momentum representation where we use the dimensionless positive wavenumber κ = kl T . The radial wavefunction, u J,n (κ), obeys the effective 1D Schrödinger equation [32,35] 4) where n = 0, 1, 2, . . . labels the solutions. For α 1, it is justified to approximate Eq. (2.4) by replacing J 2 /2κ 2 → J 2 /2α 2 . The radial problem then decouples from the angular one and becomes equivalent to a shifted 1D oscillator with energy levels (n+1/2) ω. Moreover, the angular problem reduces to a rigid rotor with the large moment of inertia α 2 / ω. We thus arrive at the E J,n quoted in Eq. (1.3), where n serves as band index and J labels the almost degenerate states within each band. We find that corrections to the energies in Eq. (1.3) scale ∼ 1/α 3 for α 1. In fact, a recent numerical study of H dot has reported that Eq. (1.3) is highly accurate for α 4 [35].
For weak-to-intermediate Coulomb interaction strength, only low-energy states are needed to span the effective Hilbert space determining the ground state. It can then be justified to retain only n = 0 modes. This step implies the restriction to angular momentum states with |J| α, since otherwise n > 0 states should also be included. For the results below, we have checked that this "single-band approximation" is indeed justified. From now on, the single-particle Hilbert space is restricted to the n = 0 sector (and the n index will be dropped). In momentum representation, this space is spanned by the orthonormal set of states [47] (2.5) Up to the zero-point contribution, the corresponding single-particle energy is E J = J 2 ω/2α 2 . The momentum-space probability density for all states is independent of J, representing a radially symmetric Gaussian peak centered at k = k 0 , The coordinate representation of Eq. (2.5) now follows by Fourier transformation, where r = r(cos θ, sin θ) with ρ = r/l T , and we use the Bessel functions J m (x) (integer m).
It will be convenient to use a second-quantized formalism below, with the noninteracting Hamiltonian which is conserved even in the interacting case. Noting that the Hamiltonian respects TRS, a finite ground-state expectation value, M s = M s = 0, corresponds to a spontaneous magnetization of the dot and thus would imply that the ground state breaks TRS. For the noninteracting case, recent work has discussed a spin-orbitinduced orbital magnetization in similar nanostructures, either in the presence [48] or absence [49] of a magnetic Zeeman field. We find below that, in the absence of a magnetic field but with strong spin-orbit coupling, already weak interactions can induce a transition to an orbital ferromagnet, where a large magnetization is present and the electrons in the dot carry a circulating charge current. This behavior appears for λ > λ c , where the critical interaction strength, λ c , vanishes in the limit α → ∞.

B. Coulomb matrix elements
The second-quantized Hamiltonian, H = H 0 + H I , with H 0 in Eq. (2.8), includes a normal-ordered Coulomb interaction term, is the Coulomb potential. Inserting the field operator (2.9), and taking into account angular momentum conservation, we find with the integer angular momentum exchange m. The real-valued Coulomb matrix elements in Eq. (2.12) take the form where we define . Using a well-known expansion formula, (2.15) where ρ > (ρ < ) is the larger (smaller) of ρ and ρ , the denominator in Eq. (2.13) is expressed as a series involving Legendre polynomials, P l (x). This allows us to perform the θ-integral in Eq. (2.13) analytically, and after some algebra we obtain with the numbers (see also Ref. [50]) and R (0) 0 = 1. The Coulomb matrix elements in Eq. (2.16) are in a convenient form for numerics [51]. In addition, as we discuss next, Eq. (2.16) also allows for analytical progress in the limit α → ∞.

C. Ultra-strong Rashba coupling
The interaction matrix elements (2.16) can be computed in closed form for α → ∞. For consistency with the single-band approximation, this limit is taken as k 0 → ∞ with l T held finite, i.e., we assume ultra-strong Rashba coupling in the presence of the dot. Taking the limit in opposite order gives similar but slightly different results; we provide a discussion of this point in the Appendix.
For α → ∞, using Eq. (2.7) and ρ = r/l T , the singleparticle states have the asymptotic real-space representatioñ (2.18) where the Gaussian e −ρ 2 /2 factor reflects the trap potential and the Rashba coupling causes rapid oscillations. Equation (2.14) is well-defined in the α → ∞ limit, G J,J+m (ρ) → π −3/2 e −ρ 2 cos(πm/2). Notably, for odd m, we find G = 0, leading to the even-odd parity effect described below. Performing the remaining integrations in Eq. (2.16), we obtain a surprisingly simple result for the Coulomb matrix elements, In terms of the R with the coefficients The small parameter η 1 in Eq. (2.20) (we take η = 0.01 for concreteness below) regularizes the l-summation, which for η = 0 is logarithmically divergent with respect to the upper limit. In physical terms, this weak divergence comes from the singular r → 0 behavior of the 1/r Coulomb potential, which in practice is cut off by the transverse (2D electron gas) confinement. Expressing the corresponding length scale as ηl T , we arrive at the regularized form in Eq. (2.20). Numerical results for the S m are shown in Table I: S m has a maximum for m = 0 and then decays with increasing |m| [52].
It is worth pointing out that the α → ∞ Coulomb matrix elements in Eq. (2.19) are valid for arbitrary radially symmetric confinement, where different confinement potentials only lead to different coefficients C l . While Eq. (2.21) describes the parabolic trap, taking for instance a hard-wall circular confinement [53], we find An ) are completely independent of the "incoming" angular momenta J 1 and J 2 . This can be rationalized by noting that in the α → ∞ limit, we arrive at an effectively homogeneous 1D problem corresponding to a ring in momentum space, see also the Appendix. For a homogeneous electron gas, on the other hand, it is well known that interaction matrix elements only depend on the exchanged (angular) momentum but not on particle momenta themselves [54]. With H 0 in Eq. (2.8), the conserved particle number, N = J c † J c J , and noting that S m = 0 for odd m, the α → ∞ Hamiltonian takes the form Since S 0 enters only via this energy shift, but otherwise disappears in H ∞ , it is convenient to put S 0 = 0 from now on and let the sum in Eq. (2.22) include m = 0; the energy E s will be kept implicit in what follows. Corrections to H ∞ at finite α originate from Coulomb matrix element contributions that vanish for α → ∞, in particular those with odd m. In Sec. III, we shall discuss the exact ground state of H ∞ for N = 2.

D. General properties of Coulomb matrix elements
We proceed by presenting symmetry relations relating different Coulomb matrix elements in Eq. (2.16). Note that our discussion here is not restricted to α → ∞, but applies to finite Rashba couplings with α 1. First, by virtue of particle indistinguishability, (2.23) Additional symmetry relations follow from the time reversal invariance of the interaction Hamiltonian H I . Indeed, because of TRS, Eq. (2.14) yields In particular, for odd m and arbitrary J, Eq. (2.24) yields The parity effect found in Sec. II B for α → ∞, with V Numerical results for α = 10 and several m are shown in Figs. 1 and 2. We draw the following conclusions: • With increasing |m|, the absolute magnitude of the Coulomb matrix elements quickly decreases.
• Pronounced differences between even and odd m are not yet visible for α = 10. Additional calculations for α = 15 and α = 30 (not shown here) confirm that the matrix elements for odd m become more and more suppressed relative to the even-m case. However, the ideal parity effect, where all odd-m matrix elements vanish for α → ∞, is approached rather slowly.
• For α = 10, Figs. 1 and 2 show that the V (m) J1,J2 carry a significant dependence on the indices (J 1 , J 2 ). This dependence ultimately disappears for α → ∞.   • For given value of m, the matrix element V (m) J1J2 has maximal absolute magnitude along the two lines J 2 = −J 1 and J 2 = J 1 + m in the (J 1 , J 2 ) plane. Noting that the single-particle eigenfunctions are localized near a ring of radius k 0 in momentum space, these two lines can be interpreted as BCSlike and exchange-type scattering processes, respec-tively, cf. the Appendix. The two lines of maximal absolute magnitude are orthogonal to each other, and cross at the point (−m/2, m/2) in the (J 1 , J 2 ) plane. While for even m, this point is not a physically realized one (since J 1,2 must be half-integer), it is always the symmetry center. , but includes the kinetic term H 0 . We assume that the interaction strength is finite, but λ 1 is needed to validate the single-band approximation. A. Two-particle eigenstates The two-particle Hilbert space is spanned by c † J1 c † J2 |0 , where we set J 1 > J 2 to avoid double counting and |0 is the N = 0 state. This space is composed of decoupled subspaces, which are invariant under the action of H ∞ . The corresponding states, |M, γ , are labeled by the integer M and a "family" index γ = 1, 2, 3, see Fig. 3 for an illustration. With amplitudes β J>0 subject to the normalization condition J>0 |β J | 2 = 1, and employing an auxiliary index i γ with values i γ=1 = 0 and i γ=2,3 = 1, those states are defined as looking for the ground-state energy, we note that an M -dependence can only originate from the E J ∼ 1/α 2 terms. For α → ∞, all |M, γ states with different M but the same γ, therefore, have the same energy. As a consequence, the interacting ground state is highly degenerate for α → ∞. This degeneracy is only lifted by finite-α corrections resulting from the kinetic energy and from Coulomb matrix elements beyond H ∞ .
Importantly, since the energy-lowering contribution −S J+J is absent in Eq. (3.2) for γ = 1, the ground state must be in the γ = 1 sector. The γ = 2, 3 states are separated by an energy gap ∼ λ ω, and we neglect these higher energy states from now on (and omit the γ index). Since the magnetization operatorM s in Eq. (2.10) is conserved, the |M states are also magnetization eigenstates. Indeed, one immediately finds that the corresponding eigenvalue is M s = 2M .

B. Distribution function
For given total angular momentum, M s = 2M , we found in Sec. III A that the eigenstate of H ∞ with lowest energy can be constructed from the Ansatz  Since the matrix appearing in Eq. (3.5) is real symmetric, we can choose real-valued β J . Moreover, since the matrix is independent of M , its lowest eigenvalue, E min , is also M -independent and depends on the interaction strength and on the Rashba coupling only through the combination λα 2 . The corresponding normalized eigenvector is easily obtained numerically and directly gives the β J . Thereby we also obtain the normalized ground-state distribution function, n J = |β J | 2 . Typical results for β J and n J are shown in Fig. 4. We find a rather broad distribution function n J , very different from a Fermi function. To reasonable approximation, the numerical results can be fitted to a Gaussian decay, n J ∼ e −(J/J * ) 2 , with J * ∼ √ α. Since J * α, the relevant angular momentum states have |J| α, and the single-band approximation is self-consistently fulfilled. As shown in the inset of Fig. 4, the β J exhibit a pairwise oscillatory behavior, where β J < 0 for J = 1/2 and 3/2, but β J > 0 for J = 5/2 and 7/2, and so on.

C. Ground state magnetization
The above results indicate that for α → ∞ and given M , the lowest energy is While this suggests that the ground state has M = 0, the M 2 /α 2 term (due to H 0 ) is in fact subleading to Coulomb corrections beyond H ∞ , which approximately scale ∼ 1/α, see Eq. (A7). We therefore have to take these Coulomb matrix elements into account when determining the ground state. To that end, using the symmetry relation (2.23), and exploiting thatM s is conserved, we note that H I [Eq. (2.12)] has the matrix elements Therefore, the energies E (∞) M in Eq. (3.6) will be independently shifted by this perturbation, and the (J 1 , J 2 )dependence of the Coulomb matrix elements becomes important, see Sec. II D. In particular, terms with odd angular momentum exchange m will contribute. Treating the Coulomb corrections in perturbation theory, the lowest energy for fixed M is  Assuming M 0 1, the main contribution to the integral comes from ϕ 1/M 0 , and performing the subsequent integration implies M 0 ≈ (λα) 1/4 . Clearly, this suggests that M 0 can be very large even for weak interactions.
The ground state of the N = 2 dot, is spanned by the two degenerate states | ± M 0 , with magnetization M s = ±2M 0 , respectively. Unless c + c − = 0, we note that |Φ is not an eigenstate of the conserved operatorM s . However, the magnetization expectation value, Φ|M s |Φ = 2(|c + | 2 −|c − | 2 )M 0 , is finite except when |c + | = |c − |. This suggests that by application of a weak magnetic field perpendicular to the 2D plane, the magnetization can be locked to one of the two minima, say, M s = +2M 0 . Adiabatically switching off the magnetic field, we then expect |Φ = |M 0 , since there is an energy barrier to the | − M 0 state. Since the barrier is not infinite, we cannot exclude that quantummechanical tunneling effects will ultimately establish an unmagnetized ground state with |c + | = |c − |, in particular when taking into account violations of the perfect rotational symmetry of the dot assumed in our model. For instance, such imperfections correspond to an eccentricity of the confinement potential or the presence of nearby impurities. As long as such imperfections represent a weak perturbation, however, the associated tunneling timescales connecting the two degenerate states | ± M 0 are expected to be very long. We shall discuss this issue in some detail in Sec. V. For practical purposes, assuming that quantum tunneling is not relevant on the timescales of experimental interest, adiabatically switching off the magnetic field then effectively results in the ground state |Φ = |M 0 . This state carries a large magnetization, M s = 2M 0 , and thus also a circulating charge current. Such a state appears to spontaneously break the TRS of the Hamiltonian, and is interpreted here as a "molecular" orbital ferromagnet.
The above discussion pertains to the idealized T = 0 case. In practice, the zero-temperature limit also governs the physics at temperatures well below the above energy barrier, k B T |δE M0 |. However, at higher temperatures, thermally induced transitions between both minima happen on short timescales, and the overall magnetization of the dot vanishes. Nonetheless,M 2 s still has a finite expectation value.

D. Spin and charge density
Before proceeding with a discussion of numerical results for N > 2, let us briefly address the spin and charge density for α → ∞. We assume that the N = 2 system is in a definite ground state, say |M 0 .
The total spin density at position r = r(cos θ, sin θ) follows as S(r, θ) = J>0 n J [s J+M0 (r, θ) + s −J+M0 (r, θ)] , (3.14) where s J = (s x J , s y J , s z J ) is the spin density for the singleparticle stateψ J (r, θ). Using Eq. (2.18), we obtain, e.g., As a consequence, the two contributions in Eq. (3.14) precisely cancel each other, and S x = 0. By the same argument, we also find that the y-and z-components of the spin density vanish. In the limit α → ∞, the spin density S is therefore identically zero. In practice, finite contributions may come from subleading (∼ 1/α) terms, but these are small for α 1. We now turn to the charge density, ρ c (r), which is always radially symmetric. For α → ∞, all single-particle states,ψ J in Eq. (2.18), lead to the same probability density in space, and we therefore conclude that ρ c (r) must be independent of Coulomb interactions. For λ 1 and arbitrary particle number N , we thus obtain ρ c (r) = eN π 3/2 l T r e −r 2 /l 2 T , (3.16) which satisfies the expected normalization, 2π ∞ 0 drrρ c (r) = eN . We mention in passing that the "edge" state property of the single-particle states, i.e., states with larger |J| live further away from the dot center, can be seen from the finite-α wavefunctions in Eq. (2.7) [32], but not anymore from their asymptotic α → ∞ form in Eq. (2.18). The λ-independence of ρ c (r) at large α is in marked contrast to the case of weak spin-orbit coupling, where ρ c contains information about interactions and can be used to detect Wigner molecule formation [5,14]. Instead, the charge density in Eq. (3.16) is featureless for arbitrary N , pointing once again to the absence of the Wigner molecule for α 1 and λ 1. Finally, we note that by computing the pair distribution function [54], we also find no trace of Wigner molecule formation in this limit.

IV. EXACT DIAGONALIZATION AND HARTREE-FOCK CALCULATIONS
We now discuss numerical results for the ground-state energy and magnetization for N ≤ 10 electrons in the dot. These results were obtained by means of the standard exact diagonalization technique and by Hartree-Fock (HF) theory from H = H 0 + H I , with H 0 in Eq. (2.8), H I in Eq. (2.12), and the Coulomb matrix elements (2.16), see Sec. II. This implies that the following results are not restricted to the α → ∞ limit considered in Sec. III. However, we are limited to rather weak interactions, λ 1, and moderate particle numbers, N < α, because of our single-band approximation. We first describe our exact diagonalization results and then turn to HF theory.

A. Exact diagonalization
Using the Rashba parameter α = 10, exact diagonalization results for N = 2 and N = 3 electrons in the dot are shown in Fig. 6. While E 0 (λ) at first sight seems rather featureless (top panel), there are non-analytic features that become visible when plotting the first derivative (center panel). Let us discuss this point in detail for N = 2, see the left side of Fig. 6. The first non-analytic feature occurs at λ c ≈ 0.25, where the second derivative diverges, d 2 E 0 /dλ 2 → −∞, as the interaction parameter λ approaches the critical value λ c from below. In close analogy to the results obtained from H ∞ in Sec. III, the ground state for λ > λ c has the magnetization M s = ±2M 0 , where M 0 is integer and the ground state is degenerate with respect to both signs. In the exact diagonalization, the "initial conditions" selecting the even- tually realized state (M s = +2M 0 or M s = −2M 0 ) correspond to unavoidable numerical rounding errors. In contrast to the α → ∞ limit, however, the interaction parameter λ must now exceed a critical value, λ c , to allow for the orbital ferromagnet. For λ < λ c , the M = 0 state is the ground state, which is adiabatically connected to the noninteracting ground state. Since energy levels of states with different conserved total angular momentum M s can cross each other, the critical value λ c marks a quantum phase transition. However, once disorder or eccentricity of the quantum dot are present, angular momentum conservation breaks down and the transition will correspond to a smooth crossover phenomenon. The observed large value of the magnetization, |M s | = 6 for λ > λ c , see Fig. 6, again rules out a purely spin-based explanation. In fact, additional jumps to even higher |M s | are observed for larger λ in Fig. 6. Similar features are also observed for N = 3, where exact diagonalization results are shown on the right side of Fig. 6. Again, the first derivative of E 0 (λ) displays non-analytic behavior. For small λ, the state stays close to a doubly degenerate Fermi sea, see Sec. II A. For λ > λ c ≈ 0.31, however, a large magnetization emerges, |M s | = 11.5 .
The results of Sec. III show that the critical interaction strength λ c vanishes for α → ∞. We therefore expect λ c to decrease with increasing α. To study this point, exact diagonalization results for α = 15 are shown in Fig. 7. All qualitative features observed for α = 10 are recovered, and the critical value λ c is indeed found to decrease: For . This confirms that with increasing spinorbit coupling strength, the orbital ferromagnetic state is reached already for weaker interactions.

B. Hartree-Fock calculations
Finally, let us turn to numerical results for larger N , where exact diagonalization becomes computationally too expensive. We have carried out an unrestricted Hartree-Fock analysis following the textbook formulation [54], in order to find the energy and the total angular momentum of the N -electron ground state. We note that (4.1) The self-consistent HF ground state is numerically found by iteration, starting from randomly chosen initial distributions. The converged {n J } distribution yields the magnetization, M s , and the ground-state energy. For λ approaching the (HF value of the) critical interaction parameter, λ c (α, N ), the energy shows similar nonanalytic features as found from exact diagonalization, see Sec. IV A. For λ > λ c , a large ground-state magnetization is observed, again corresponding to orbital ferromagnetism.
Our HF results for λ c and M s are shown in Fig. 8. We consider the Rashba spin-orbit coupling α = 30, and up to N = 10 electrons in the dot. Unfortunately, we cannot address larger N , for otherwise our single-band approximation is not justified anymore. For N = 2, the corresponding exact diagonalization values are also given. The HF value for λ c is only slightly smaller than the exact one, which suggests that HF theory is at least qualitatively useful. That the HF prediction is below the exact one for N = 2 can be rationalized by noting that HF theory generally tends to favor ordered phases such as orbital ferromagnetism, resulting in a smaller value for λ c . The magnetization for λ > λ c , however, is a more difficult quantity to predict due to the shallow minima of the free energy curves in Fig. 5. Indeed, the inset of Fig. 8 shows that the HF value of the magnetization (which appears to scale as M s ∼ N ) is significantly smaller than the exact one. With increasing N , the HF predictions for λ c indicate that the transition to the orbital ferromagnet persists. Moreover, this transition can even be reached at weaker interactions.

V. DISCUSSION
In this work, we have studied the interacting Nelectron problem for a parabolic 2D quantum dot (with N ≤ 10) in the limit of strong Rashba spin-orbit coupling, α 1. This regime is characterized by an almost flat single-particle spectrum, where we find that already weak-to-intermediate Coulomb interactions (our singleband approximation permits us to study the regime λ 1 only) are sufficient to induce molecular orbital ferromagnetism. This state is observed for λ > λ c (α, N ), where our N = 2 solution in Sec. II shows that the critical strength λ c → 0 for α → ∞. For finite (but large) α, however, λ c will be finite. The orbital ferromagnet has a giant total angular momentum, accompanied by a circulating charge current.
Coming back to our discussion in Sec. III C, we now address issues concerning the experimental observation of the predicted orbital ferromagnetism for a single quantum dot. The transition to this state could be induced in practice by varying the electrostatic confinement potential and/or the gate-controlled Rashba spin-orbit coupling in order to reach the regime defined by α 1 and λ > λ c . By allowing for an eccentricity of the dot confinement potential, which also can be achieved with appropriate gate voltages, quantum tunneling processes connecting the free energy minima with opposite magnetization, M s = ±M min , are expected to become relevant, see Sec. III C. The corresponding timescale for such processes can be estimated as follows. We first note that the free energy barrier between both minima, B ω, corresponds to a number B ≈ 0.1 . . . 0.15, see Fig. 5. We next employ a paradigmatic effective low-energy model to include the effects of imperfections breaking the ideal rotational symmetry, (5.1) The first term describes the dot eccentricity, with a small dimensionless parameter , where the polar angle φ is conjugate to the magnetization M s . The second term approximates the double-well potential in Fig. 5. The two lowest eigenenergies for H eff are known exactly [58]. From the result, we find the level splitting The resulting timescale for tunneling processes, τ , is thereby estimated as ωτ = ω δE ≈ 0.2e 5.96Mmin . (5.3) where, for simplicity, we have put B = 0.1 and = 0.01. For the value M min ≈ 18 observed in Fig. 5, we get the estimate ωτ ≈ 10 45 . This astronomically long tunneling time strongly suggests that on experimentally accessible timescales, the orbital ferromagnet described in this paper will be indistinguishable from a true equilibrium state. It is also useful to contrast the behavior reported here to the well-known persistent currents in normal-metal quantum rings [59][60][61][62][63], where a circulating equilibrium electric current flows and can be experimentally detected, see Ref. [64] and references therein. First, a persistent current flows already in noninteracting quantum rings but requires a nonzero flux threading the ring, while orbital ferromagnetism in a 2D dot is generated by the interplay of Coulomb interactions and strong spin-orbit coupling. Second, the total angular momentum (magnetization) predicted here for a 2D Rashba dot can be very large. Therefore, the circulating currents in our case should exceed the persistent currents observed in quantum rings by far. Despite of these differences, the persistent current analogy also suggests ways to observe our predictions experimentally. Another possibility is to study the response to a weak magnetic field applied perpendicular to the 2D plane. The low-field susceptibility is then expected to be singular, just as in an ordinary ferromagnet. At elevated temperatures approaching the free energy barrier height discussed above, the orbital magnetization in our system will be thermally suppressed and ultimately disappear. The relevant temperature scale for this crossover is T c ≈ B ω/k B . For typical quantum dots [2], T c ≈ 1 to 10 K.
To conclude, we hope that our prediction of orbital ferromagnetism in Rashba dots will stimulate further theoretical and experimental work. For instance, it remains an open question to address the transition from the orbital ferromagnet to a Wigner molecule with increasing interaction strength for large Rashba coupling. In order to achieve this description, one needs to go beyond the single-band approximation employed in this work. × δ(k 1 − k 0 ) δ(k 2 − k 0 )e im(φ2−φ1) σ1,σ2=± where the δ-function implies the constraint k 1 (φ 1 ) = k 2 0 + q 2 + 2k 0 q cos φ 1 = k 0 , and similarly for k 2 . This leads to the condition cos φ 1 = − cos φ 2 = −q/2k 0 , which is met by two types of scattering processes only, namely (a) for φ 2 = π + φ 1 (BCS-like pairing), and (b) for φ 2 = π − φ 1 (exchange-type process), see Fig. 9. Such spin-orbit-induced constraints on interaction processes were also recently pointed out in Ref. [66]. where the first (second) term in the numerator results from BCS-like (exchange-type) processes. Importantly, the above integral is infrared divergent for q = 2k 0 cos ϕ → 0. To regularize this singularity, we employ l T as effective system size and require ql T > 1. After some algebra, we find the (J 1 , J 2 )-independent result V (m) J1,J2 λ ωδ m,even , which recovers the parity effect in Sec. II C, including the (J 1 , J 2 )-independence of the matrix elements. In contrast to Eq. (2.19), however, the even-m Coulomb matrix elements found here are also independent of m. This indicates that the limits k 0 → ∞ and l T → ∞ do not commute.