Electron and boson clusters in confined geometries: Symmetry breaking in quantum dots and harmonic traps
-
Edited by A. Welford Castleman, Jr., Pennsylvania State University, University Park, PA, and approved January 26, 2006 (received for review October 20, 2005)
Abstract
We discuss the formation of crystalline electron clusters in semiconductor quantum dots and of crystalline patterns of neutral bosons in harmonic traps. In a first example, we use calculations for two electrons in an elliptic quantum dot to show that the electrons can localize and form a molecular dimer. The calculated singlet–triplet splitting (J) as a function of the magnetic field (B) agrees with cotunneling measurements with its behavior reflecting the effective dissociation of the dimer for large B. Knowledge of the dot shape and of J(B) allows determination of the degree of entanglement. In a second example, we study strongly repelling neutral bosons in two-dimensional harmonic traps. Going beyond the Gross–Pitaevskii (GP) mean-field approximation, we show that bosons can localize and form polygonal-ring-like crystalline patterns. The total energy of the crystalline phase saturates in contrast to the GP solution, and its spatial extent becomes smaller than that of the GP condensate.
Explorations of the size-dependent evolution of the properties of materials are at the frontier of modern condensed matter and materials research. Indeed, investigations of clusters containing a finite, well defined number of elementary building units (atoms, molecules, electrons, or other elementary constituents) allow investigations of the transition from the atomic or molecular regime to the finite nano-aggregate domain, and ultimately of the convergence with increasing size to the condensed phase, extended system category. Moreover, investigations of clusters provide opportunities for discovery of novel properties and phenomena that are intrinsic properties of finite systems, distinguishing them from bulk materials (1).
Commonly, studies of materials clusters involve atoms or molecules interacting through electrostatic or electromagnetic potentials, with the heavier nuclei being the “structural skeleton” and the much lighter electrons serving as the “glue” that binds the atoms together. In this article we focus on novel, somewhat exotic, types of clusters. In particular, we discuss clusters of electrons in man-made (artificial) quantum dots (QDs) created through lithographic and gate-voltage techniques at semiconductor interfaces, and clusters of neutral atoms in traps under conditions that may relate to formation of Bose–Einstein condensates (BECs). We illustrate that these cluster systems reveal interesting emergent physical behavior arising from spontaneous breaking of spatial symmetries; symmetry breaking is defined as a circumstance where a lower energy solution of the Schrödinger equation is found that is characterized by a lower symmetry than that of the Hamiltonian of the system. Such symmetry breaking is exhibited through the formation of clusters of localized electrons (often called Wigner molecules, or WMs) in 2D QDs [see Fig. 1 and in particular our discussion of two-electron WMs below (2)]. Symmetry breaking is also manifested in the transition (3), induced by increasing the interatomic repulsive interaction strength, of the BEC state of neutral atoms confined by a parabolic 2D trap to a crystalline cluster state.
UHF electron density in a parabolic QD for N = 19 electrons and S z = 19/2, exhibiting breaking of the circular symmetry at R W = 5 and zero magnetic field. Remaining parameters: parabolic confinement, ℏω0 = 5 meV; effective mass, m* = 0.067me. Distances are in nanometers, and the electron density is in 10−4 nm−2.
Two-Electron WMs
Electron localization leading to formation of molecular-like structures (the aforementioned WMs) within a single circular 2D QD at zero magnetic field (B) has been theoretically predicted to occur (4–11), as the strength of the e–e repulsive interaction relative to the zero-point energy increases, as expressed through an increasing value of the Wigner
parameter R
W, defined as R
W = Z
2
e
2/(ℏω0
l
0), with l
0 =
being the characteristic harmonic-oscillator length of the confining potential. [The subscript W in the case of a Coulomb
force stands for “Wigner,” because the confined clusters of localized electrons may be viewed as finite-size precursors of
the bulk Wigner crystal (12).] The formation of such “molecular structures” is a manifestation of spontaneous symmetry breaking associated with a quantum
phase transition, occurring at zero magnetic field for circular 2D QDs for R
W ≥ 1 and involving the crossover from a liquid-like state to a crystalline one. This crossover, described by using symmetry
breaking (4), was confirmed in other studies (5, 7–9) using a variety of methods. The degree of electron localization, which underlies the appearance of crystalline patterns,
has been described (4) as a progression from “weak” to “strong” WMs as a function of increasing R
W (or equivalently decreasing density). In high magnetic fields, the rotating electron molecules exhibit magic angular momenta
(13–16) corresponding to fractional quantum Hall effect fillings. This has led (14–16) to the derivation of an analytic trial wave function that provides a better description of the finite-size analogs of the
fractional quantum Hall effect in comparison with the Laughlin (17) and composite-fermion (18) wave functions.
Because of the finite size, WMs are expected to show new behavior that differs from the classical Wigner crystal familiar from solid state physics. The limit of a classical Wigner crystal is expected to be reached for a higher number of electrons N and very large R W (4). In the following we use the term “Wigner molecule” even in the case of only two localized electrons. For this case the WM exhibits close analogies to an H2 natural molecule, as described below.
Here, we focus on a two-electron (2e) WM, in light of the current experimental effort (19, 20) aiming at implementation of a spin-based (21) solid-state quantum logic gate that employs two coupled one-electron QDs (double dot). We present an exact diagonalization (EXD) and an approximate (generalized Heitler–London, GHL) microscopic treatment for two electrons in a single elliptic QD specified by the parameters of a recently investigated experimental device (22). While formation of WMs in circular QDs requires weak confinement (that is, small ω0 in the expression for R W given above), and thus large dots of lower densities (so that the interelectron repulsion dominates), we show that formation of such WMs is markedly enhanced in highly deformed (e.g., elliptic) dots due to their lower symmetry. The calculations provide a good description of the measured J(B) curve [the singlet–triplet (ST) splitting] when screening (23, 24) due to the metal gates and leads is included (in addition to the dielectric constant of the semiconductor, GaAs). In particular, our results reproduce the salient experimental findings pertaining to the vanishing of J(B) for a finite value of B ∼ 1.3 T [associated with a change in sign of J(B) indicating a ST transition], as well as the flattening of the J(B) curve after the ST crossing. These properties, and in particular the latter one, are related directly to the formation of an electron molecular dimer and its effective dissociation for large magnetic fields. The effective dissociation of the electron dimer is most naturally described through the GHL approximation, and it is fully supported by the more accurate, but physically less transparent, EXD.
Of special interest for quantum computing (21) is the degree of entanglement exhibited by the two-electron molecule in its singlet state. Entanglement is a purely quantum mechanical phenomenon in which the quantum state of two or more objects cannot be described independently of each other, even when the individual objects are spatially separated. The highest degree of entanglement occurs at full separation, as discussed in the celebrated Einstein, Podolsky, and Rosen paper (25). Electrons confined in a QD are not necessarily spatially separated from each other, and consequently their degree of entanglement may be lower than the maximal one, as shown by us below.
Here, in relation to the microscopic calculations, we investigate two different measures of entanglement. The first, known as the concurrence (𝒞) for two indistinguishable fermions (26, 27), has been used in the analysis of the experiment in ref. 22 [this measure is related to the operational cycle of a two-spin-qubit quantum logic gate (26, 27)]. The second measure, referred to as the von Neumann entropy (𝒮) for indistinguishable particles, was developed in ref. 28 and used in ref. 29. We show that the present wave-function-based methods, in conjunction with the knowledge of the dot shape and the J(B) curve, enable theoretical determination of the degree of entanglement, in particular for the elliptic QD of ref. 22. The increase in the degree of entanglement (for both measures) with stronger magnetic fields correlates with the dissociation of the 2e molecule. This supports the experimental assertion (22) that cotunneling spectroscopy can probe properties of the electronic wave function of the QD, and not merely its low-energy spectrum. Our methodology can be straightforwardly applied to other cases of strongly interacting devices, e.g., double dots with strong interdot tunneling.
Clusters of Neutral Bosons in Harmonic Traps
BECs in harmonic traps (30, 31) are normally associated with weakly interacting neutral atoms, and their physics is described adequately by the Gross–Pitaevskii (GP) mean-field theory (32). Lately, however, experimental advances in controlling the interaction strength (33–36) permit the production of novel bosonic states in the regime of strong interparticle repulsions. Theoretical efforts motivated by this capability include studies of the Bose–Hubbard model (37, 38), and investigations about the “fermionization” limit of a 1D gas of trapped impenetrable bosons (39–41), often referred to as the Tonks–Girardeau (TG) regime (39, 42). Here, we address the problem of strongly repelling (impenetrable) bosons in higher dimensions. In particular, we discuss 2D interacting bosons in a circular harmonic trap, with the extension to 3D systems being straightforward. To this end, we use computational methods that go beyond the GP method.
We explore the transition from a BEC (diffuse cloud) state to a crystalline phase, in which the trapped localized bosons form crystalline patterns. At the mean-field level, these crystallites are static and are portrayed directly in the single-particle densities. After restoration of rotational symmetry, the single-particle densities are circularly symmetric, and thus the crystalline symmetry becomes “hidden”; however, it can be revealed in the conditional probability distribution (CPD) (anisotropic pair correlation), P(r, r 0), which expresses the probability of finding a particle at r given, that the “observer” (i.e., reference point) is riding on another particle at r 0 (10, 16).
Methods
Two-Electron QD: Microscopic Treatment.
The Hamiltonian for two 2D interacting electrons is
where the last term is the Coulomb repulsion, κ is the dielectric constant, and r
12 = |r
1 − r
2|. H(r) is the single-particle Hamiltonian for an electron in an external perpendicular magnetic field B and an appropriate confinement potential. When position-dependent screening is included, the last term in Eq. 1 is modified by a function of r
12 (see below). For an elliptic QD, the single-particle Hamiltonian is written as
where T = (p − e
A/c)2/2m*, with A = 0.5(−By, Bx, 0) being the vector potential in the symmetric gauge. m* is the effective mass, and p is the linear momentum of the electron. The second term is the external confining potential; the last term is the Zeeman
interaction with g* being the effective g factor, μB is the Bohr magneton, and s is the spin of an individual electron.
The GHL method for solving the Hamiltonian (1) consists of two steps. In the first step, we solve selfconsistently the ensuing unrestricted Hartree–Fock (UHF) equations allowing for lifting of the double-occupancy requirement (imposing this requirement gives the restricted HF method, RHF). For the S z = 0 solution, this step produces two single-electron orbitals u L,R(r) that are localized left (L) and right (R) of the center of the QD [unlike the RHF method that gives a single doubly-occupied elliptic (and symmetric about the origin) orbital]. At this step, the many-body wave function is a single Slater determinant ΨUHF (1 ↑, 2 ↓) ≡ | u L(1 ↑)u R(2 ↓)〉 made out of the two occupied UHF spin-orbitals uL(1 ↑) ≡ uL(r 1)α (1) and u R(2 ↓) ≡ u R(r 2)β(2), where α(β) denotes the up (down) [↑ (↓)] spin. This UHF determinant is an eigenfunction of the projection S z of the total spin S = s 1 + s 2 but not of S 2 (or the parity space-reflection operator).
In the second step, we restore the broken parity and total-spin symmetries by applying to the UHF determinant the projection
operator (43, 44) P
s
,t = 1 ∓ ϖ12, where the operator ϖ12 interchanges the spins of the two electrons; the upper (minus) sign corresponds to the singlet. The final result is a GHL
two-electron wave function ΨGHL
s,t(r
1, r
2) for the ground-state singlet (index s) and first-excited triplet (index t), which uses the UHF localized orbitals,
where χs
,t = (α(1)β(2) ∓ α(2)β(1)) is the spin function for the 2e singlet and triplet states. The general formalism of the 2D UHF equations and of the subsequent restoration of broken spin
symmetries can be found in refs. 11 and 43–45.
The use of optimized UHF orbitals in the GHL is suitable for treating single elongated QDs. The GHL is equally applicable to double QDs with arbitrary interdot-tunneling coupling (43, 44). In contrast, the Heitler–London treatment (46) (known also as valence bond), where nonoptimized “atomic” orbitals of two isolated QDs are used, is appropriate only for the case of a double dot with small interdot-tunneling coupling (21).
The orbitals u
L,R(r) are expanded in a real Cartesian harmonic-oscillator basis, i.e.,
where the index j ≡ (m, n) and ϕj(r) = Xm(x)Yn(y), with Xm(Yn) being the eigenfunctions of the 1D oscillator in the x(y) direction with frequency ωx(ωy). The parity operator 𝒫 yields 𝒫Xm(x) = (−1)mXm(x), and similarly for Yn(y). The expansion coefficients C
j
L,R are real for B = 0 and complex for finite B. In the calculations, we use K = 79, yielding convergent results.
In the exact-diagonalization method, the many-body wave function is written as a linear superposition over the basis of noninteracting
two-electron determinants, i.e.,
where ψ(1; i) = ϕi(1 ↑) if 1 ≤ i ≤ K and ψ(1; i) = ϕi
−K(1 ↓) if K + 1 ≤ i ≤ 2K [and similarly for ψ(2; j)]. The total energies E
j
L,R and the coefficients Ωij
s,t are obtained through a “brute force” diagonalization of the matrix eigenvalue equation corresponding to the Hamiltonian in
Eq. 1. The EXD wave function does not immediately reveal any particular form, although our calculations below show that it can
be approximated by a GHL wave function in the case of the elliptic dot under consideration.
Two-Electron QD: Measures of Entanglement.
To calculate the concurrence 𝒞 (26, 27), one needs a decomposition of the GHL wave function into a linear superposition of orthogonal Slater determinants. Thus, one needs to expand the nonorthogonal
u
L,R(r) orbitals as a superposition of two other orthogonal ones. To this effect, we write u
L,R(r) ∝ Φ+(r) ± ξΦ−(r), where Φ+(r) and Φ−(r) are the parity symmetric and antisymmetric (along the x axis) components in their expansion given by Eq. 4. Subsequently, with the use of Eq. 3, the GHL singlet can be rearranged as follows:
where the so-called interaction parameter (27), η = ξ2, is the coefficient in front of the second determinant. Knowing η allows a direct evaluation of the concurrence of the singlet
state, since 𝒞s = 2η/(1 + η2) (27). Note that Φ+(r) and Φ−(r) are properly normalized. It is straightforward to show that η = (1 − |S
LR|)/(1 + |S
LR|), where S
LR (with |S
LR| ≤ 1) is the overlap of the original u
L,R(r) orbitals.
For the GHL triplet, one obtains an expression independent of the interaction parameter η, i.e.,
which is a maximally (𝒞t = 1) entangled state. Note that underlying the analysis of the experiments in ref. 22 is a conjecture that wave functions of the form given in Eqs. 6 and 7 describe the two electrons in the elliptic QD.
To compute the von Neumann entropy, one needs to bring both the EXD and the GHL wave functions into a diagonal form [the so-called
“canonical form” (28, 47)], i.e.,
with the Φ(i)s being appropriate spin orbitals resulting from a unitary transformation of the basis spin orbitals ψ(j)s (see Eq. 5); only terms with zk ≠ 0 contribute. The upper bound M can be smaller (but not larger) than K (the dimension of the single-particle basis); M is referred to as the Slater rank. One obtains the coefficients of the canonical expansion from the fact that the |zk|2 are eigenvalues of the hermitian matrix Ω†Ω [Ω (see Eq. 5) is antisymmetric]. The von Neumann entropy is given by
with the normalization Σk=1
M |zk|2 = 1. Note that the GHL wave functions in Eqs. 6 and 7 are already in canonical form, which shows that they always have a Slater rank of M = 2. One finds 𝒮GHL
s = log2(1 + η2) − η2 log2(η2)/(1 + η2), and 𝒮GHL
t = 1 for all B. For large B, the overlap between the two electrons of the dissociated dimer vanishes, and thus η → 1 and 𝒮GHL
s → 1.
Neutral Repelling Bosons in Harmonic Traps: GP Equation.
Mean-field symmetry breaking for bosonic systems has been discussed earlier in the context of two-component condensates, where
each species is associated with a different space orbital (48, 49). We consider here one species of bosons but allow each particle to occupy a different space orbital φi(r
i). The permanent |ΦN〉 = Perm[φ1(r
1), …, φN(r
N)] serves as the many-body wave function of the unrestricted Bose–Hartree–Fock (UBHF) approximation. This wave function reduces to the GP form with the restriction that all bosons occupy the same orbital φ0(r), i.e., |ΦN
GP〉 = Πi=1
Nφ0(r
i), and φ0(r) is determined self-consistently at the restricted Bose–Hartree–Fock (RBHF) level via the equation (50)
Here, V(r
1, r
2) is the two-body repulsive interaction, which is taken to be a contact potential, V
δ = gδ(r
1 − r
2), for neutral bosons. The single-particle Hamiltonian is given by H
0(r) = − ℏ2▿2/(2m) + mω0
2
r
2/2, where ω0 characterizes the harmonic confinement.
Neutral Repelling Bosons in Harmonic Traps: Symmetry Breaking.
We simplify the solution of the UBHF problem by considering explicit analytic expressions for the space orbitals φi(r i). In particular, since the bosons must avoid occupying the same position in space to minimize their mutual repulsion, we take all of the orbitals to be of the form of displaced Gaussians, namely, φi(r i) = π−1/2σ−1 exp[−(r i − a i)2/(2σ2)]. The positions a i describe the vertices of concentric regular polygons, with both the width σ and the radius a = |a i| of the regular polygons determined variationally through minimization of the total energy E UBHF = 〈ΦN|H|ΦN〉/〈ΦN|ΦN〉, where H = Σi=1 N H 0(r i) + Σi<j N V(r i, r j) is the many-body Hamiltonian.
With the above choice of localized orbitals, the unrestricted permanent |ΦN〉 breaks the continuous rotational symmetry. However, the resulting energy gain becomes substantial for stronger repulsion. Controlling this energy gain (the strength of correlations) is the ratio R δ between the strength of the repulsive potential and the zero-point kinetic energy. Specifically, for a 2D trap, one has R δ = gm/(2π ℏ2) for a contact potential.
Neutral Repelling Bosons in Harmonic Traps: Restoration of Broken Symmetry.
Although the optimized UBHF permanent |ΦN〉 performs exceptionally well regarding the total energies of the trapped bosons, in particular in comparison with the restricted
wave functions (e.g., the GP anzatz), it is still incomplete. Indeed, due to its localized orbitals, |ΦN〉 does not preserve the circular (rotational) symmetry of the 2D many-body Hamiltonian H. Instead, it exhibits a lower point-group symmetry, i.e., a C
2 symmetry for N = 2 and a C
5 one for the (1, 5) structure of N = 6 (see below). As a result, |ΦN〉 does not have a good total angular momentum. This paradox is resolved through a post-Hartree–Fock step of restoration of broken symmetries via projection techniques (45, 51), yielding a new wave function |ΨN,L
PRJ〉 with a definite angular momentum L, that is
where |ΦN(γ)〉 is the original UBHF permanent having each localized orbital rotated by an azimuthal angle γ, with L being the total angular momentum. The projection yields wave functions for a whole rotational band. Note that the projected
wave function |ΨN,L
PRJ〉 in Eq. 11 may be regarded as a superposition of the rotated permanents |ΦN(γ)〉, thus corresponding to a “continuous-configuration-interaction” solution.
Here, we are interested in the projected ground-state (L = 0) energy, which is given by
Results and Discussion
Two-Electron QD.
To model the experimental elliptic QD device, we take, following ref. 22, ℏωx = 1.2 meV and ℏωy = 3.3 meV. The effective mass of the electron is taken as m* = 0.067 me (GaAs). Since the experiment did not resolve the lifting of the triplet degeneracy caused by the Zeeman term, we take g* = 0. Using the two-step method, we calculate the GHL ST splitting J GHL(B) = E GHL s(B) − E GHL t(B) as a function of the magnetic field in the range 0 ≤ B ≤ 2.5 T. Screening of the e–e interaction due to the metal gates and leads must be considered to reproduce the experimental J(B) curve. This screening can be modeled, to first approximation, by a position-independent adjustment of the dielectric constant κ (52). Indeed, with κ = 22.0 (instead of the GaAs dielectric constant; i.e., κ = 12.9), good agreement with the experimental data is obtained (see Fig. 2). In particular, we note the ST crossing for B ≈ 1.3 T and the flattening of the J(B) curve beyond this crossing.
The ST splitting J = Es − Et as a function of the magnetic field B for an elliptic QD with ℏωx = 1.2 meV and ℏωy = 3.3 meV (these values correspond to the device of ref. 22). Solid line, GHL (broken-symmetry UHF plus restoration of symmetries) results with a coordinate-independent screening (κ = 22); dashed line, EXD results with κ = 12.9 (GaAs) but including screening with a coordinate dependence according to ref. 24 and d = 18.0 nm (see text). Remaining material parameters: m*(GaAs) = 0.067m e, and g* = 0 (see text). The experimental measurements (22) are denoted by open squares. Our sign convention for J is opposite of that in ref. 22.
We have also explored, particularly in the context of the EXD treatment, a position-dependent screening using the functional form, (e 2/κr 12)[1 − (1 + 4d 2/r 12 2)−1/2], proposed in ref. 24, with d as a fitting parameter. The J EXD(B) result for d = 18.0 nm is depicted in Fig. 2 (dotted line) and is in very good agreement with the experimental measurement.
The singlet-state electron densities from the GHL and the EXD treatments at B = 0 and B = 2.5 T are displayed in Fig. 3. These densities illustrate the dissociation of the electron dimer with increasing magnetic field. The asymptotic convergence (beyond the ST point) of the energies of the singlet and triplet states, i.e., [J(B) → 0 as B → ∞], is a reflection of the dissociation of the 2e molecule, since the ground-state energy of two fully spatially separated electrons (zero overlap) does not depend on the total spin.
Total electron densities (EDs) associated with the singlet state of the elliptic dot at B = 0 and B = 2.5 T. (a) The GHL densities. (b) The EXD densities. The rest of the parameters and the screening of the Coulomb interaction are as in Fig. 2. Lengths are in nanometers, and densities are in 10−4 nm−2.
In contrast, the singlet-state RHF electron densities fail to exhibit formation of an electron dimer for all values of B. This underlies the failure of the RHF method to describe the behavior of the experimental J(B) curve. In particular, J RHF(B = 0) has the wrong sign, while J RHF(B) diverges for high B as is the case for the RHF treatment of double dots (see ref. 43).
For the GHL singlet, using the overlaps of the left and right orbitals, we find that starting with η = 0.46 (𝒞 s = 0.76) at B = 0, the interaction parameter (singlet-state concurrence) increases monotonically to η = 0.65 (𝒞 s = 0.92) at B = 2.5 T. At the intermediate value corresponding to the ST transition (B = 1.3 T), we find η = 0.54 (𝒞 s = 0.83).
Our B = 0 theoretical results for η and 𝒞 s are in remarkable agreement with the experimental estimates (22) of η = 0.5 ± 0.1 and 𝒞 s = 0.8, which were based solely on conductance measurements below the ST transition (i.e., near B = 0). We note that, for the RHF, 𝒞 RHF s = 0, since a single determinant is unentangled for both the two measures considered here.
Since the EXD singlet has obviously a Slater rank M > 2, the definition of concurrence is not applicable to it. The von Neumann entropy for the EXD singlet (𝒮EXD s) is displayed in Fig. 4, along with that (𝒮GHL s) of the GHL singlet. 𝒮EXD s and 𝒮GHL s are rather close to each other for the entire B range, and it is remarkable that both remain close to unity for large B, although the maximum allowed mathematical value is log2(K) [as aforementioned, we use K = 79; i.e., log2(79) = 6.3]; this maximal value applies for both the EXD and GHL approaches. The saturation of the entropy for large B to a value close to unity reflects the dominant (and roughly equal at large B) weight of two configurations in the canonical expansion (see Eq. 8) of the EXD wave function, which are related to the two terms (M = 2) in the canonical expansion of the GHL singlet (Eq. 6). This is illustrated by the histograms of the |z k s|2 coefficients for B = 1.3 T at the top of Fig. 4. These observations support the GHL approximation, which is computationally less demanding than the exact diagonalization, and can be used easily for larger N.
Von Neumann entropy for the singlet state of the elliptic dot as a function of the magnetic field B. Solid line, GHL; dashed line, EXD. The rest of the parameters and the screening of the Coulomb interaction are as in Fig. 2. At the top, we show histograms for the |zk|2 coefficients (see Eq. 8) of the singlet state at B = 1.3 T, illustrating the dominance of two configurations. Note the small third coefficient |z 3|2 = 0.023 in the EXD case.
Neutral Repelling Bosons in Harmonic Traps.
In Fig. 5, we display as a function of the parameters R δ the total energies for N = 6 bosons calculated at several levels of approximation. In both cases the lowest UBHF energies correspond to a (1,5) crystalline configuration, namely one boson is at the center and the rest form a regular pentagon of radius a. Observe that the GP total energies are slightly lower than the E RBHF G ones; however, both exhibit an unphysical behavior since they diverge as R δ → ∞. This behavior contrasts sharply with that of the UHF energies, E UBHF and PRJ (see below), which saturate as R δ → ∞; in fact, a value close to saturation is achieved already for R δ ∼ 10. We have checked that for all cases with N = 2 − 7, the total energies exhibit a similar behavior. For a repulsive contact potential, the saturation of the UBHF energies is associated with the ability of the trapped bosons (independent of N) to minimize their mutual repulsion by occupying different positions in space, and this is one of our central results. For N = 2, the two bosons localize at a distance 2a apart to form an antipodal dimer. For N ≤ 5, the preferred UBHF crystalline arrangement is a single ring with no boson at the center [usually denoted as (0, N)]. N = 6 is the first case having one boson at the center [designated as (1, N − 1)], and the (0, 6) arrangement is a higher energy isomer. The structural parameters (e.g., the width of the Gaussian orbitals and the radii of the polygonal ring, calculated via the UBHF method) show a saturation behavior similar to that illustrated above for the energy of the system (3). In contrast, the width of the condensate cloud (i.e., the GP solution) diverges with increasing repulsion strength (R δ).
Total energies as a function of R δ for various approximation levels, calculated for N = 6 harmonically confined 2D bosons in the (1,5) lowest-energy configuration. RBHF/G, restricted Bose–Hartree–Fock (RBHF) energy, E RBHF G , with the common orbital φ0(r) approximated by a Gaussian centered at the trap origin; GP, the Gross–Pitaevskii energy; PRJ, the energy of the symmetry-restored state obtained via projection of the (unrestricted) UBHF state. Energies are in units of ℏω0.
The saturation found here for 2D trapped bosons interacting through strong repelling contact potentials is an illustration of the “fermionization” analogies that appear in strongly correlated systems in all three dimensionalities. Indeed such energy saturation has been shown for the TG 1D gas (42, 39) and has also been discussed for certain 3D systems [i.e., three trapped bosons (53) and an infinite boson gas (54)]. Saturation of the energy and the length of the trapped atom cloud (and thus of the interparticle distance) has been measured recently for the 1D TG gas (see in particular figures 3 and 4 in ref. 36, and compare with the similar trends predicted here for the 2D case in Fig. 5).
For N = 6 2D bosons, Fig. 5 shows that the E 0 PRJ energies share with the UBHF ones the saturation property for the case of a contact-potential repulsion. However, the projection brings further lowering of the total energies compared with the UBHF ones. [The projected ground state is always lower in energy than the original broken-symmetry one (55).] Thus, for strong interactions (large values of R δ) the restoration-of-broken-symmetry step yields an excellent approximation of both the exact many-body wave function and the exact total energy.
The transformations of the single-particle densities (displayed in Fig. 6 for N = 6 neutral bosons interacting via a contact potential and R δ = 25) obtained from application of the successive approximations provide an illustration of the two-step method of symmetry breaking with subsequent symmetry restoration. Indeed, the GP single-particle density (Fig. 6 a) is circularly symmetric, but the UBHF one (Fig. 6 b) explicitly exhibits a (1,5) crystalline configuration. After symmetry restoration (Fig. 6 c) the circular symmetry is reestablished, but the single-particle density is radially modulated unlike the GP density. In addition, the crystalline structure in the projected wave function is now hidden; however, it can be revealed through the use of the conditional probability distribution (10, 16) (see Fig. 6 d), which resembles the (crystalline) UBHF single-particle density, but with one of the humps on the outer ring missing (where the observer is located). In particular, P(r 0, r 0) ≈ 0 and the boson associated with the observer is surrounded by a “hole” similar to the exchange-correlation hole in electronic systems. This is another manifestation of the “fermionization” of the strongly repelling 2D bosons. However, here as in the 1D TG case (42, 39), the vanishing of P(r 0, r 0) results from the impenetrability of the bosons. For the GP condensate, the conditional probability distribution is independent of r 0, i.e., P GP(r, r 0) ∝ |φ0(r)|2, reflecting the absence of any space correlations.
Single-particle densities (a–c) and conditional probability distribution (d) for N = 6 2D harmonically trapped neutral bosons with a contact interaction and R δ = 25. (a) The single-orbital self-consistent GP case. (b) The symmetry-broken UBHF case (static crystallite). (c) The projected case (symmetry-restored wave function; see Eq. 11). The crystalline structure of the outer ring in this last case is “hidden” but is revealed in the conditional probability distribution (10, 16) displayed in d, where the observation point is denoted by a black dot (on the right). Lengths are in units of l 0.
It is of importance to observe that the radius of the BEC (GP case; Fig. 6 a) is significantly larger than the actual radius of the strongly-interacting crystalline phase (projected wave function; Fig. 6 c). This is because the extent of the crystalline phase saturates, whereas that of the GP condensate grows with no bounds as R δ → ∞. Such dissimilarity in size (between the condensate and the strongly interacting phase) has been also predicted (40) for the trapped 1D Tonks–Girardeau gas and indeed observed experimentally (36). In addition, the 2D single-particle momentum distributions for neutral bosons have a one-hump shape with a maximum at the origin (a behavior exhibited also by the trapped 1D TG gas). The width of these momentum distributions versus R δ increases and saturates to a finite value, whereas that of the GP solution vanishes as R δ → ∞.
Summary
In this article we explored symmetry-breaking transitions predicted to occur in confined fermionic and bosonic systems when the strength of the interparticle repulsive interactions exceeds an energy scale that characterizes the degree of confinement. For two electrons in an elliptic QD, we predicted formation and effective dissociation (with increasing magnetic field) of an electron dimer, which is reflected in the behavior of the computed ST splitting, J(B), that agrees well (Fig. 2) with measurements (22).
Furthermore, we showed that, from a knowledge of the dot shape and of J(B), theoretical analysis along the lines introduced here allows probing of the correlated ground-state wave function and determination of its degree of entanglement. This presents an alternative to the experimental study where determination of the concurrence used conductance data (22). Such information is of interest to the implementation of spin-based solid-state quantum logic gates.
For the case of 2D trapped bosonic clusters, we found with increasing repulsive two-body interaction localization of the bosons in the trap, resulting in formation of crystalline patterns made of polygonal rings; although we have focused here on repulsive contact interactions, similar results were obtained also for a Coulomb repulsion (3).
These results provide the impetus for experimental efforts to access the regime of strongly repelling neutral bosons in two dimensions. To this end we anticipate that extensions of methodologies developed for the recent realization of the Tonks–Girardeau regime in one dimension [using a finite small number of trapped 87Rb and optical lattices, with a demonstrated wide variation of R δ from 5 to 200 (35) and from 1 to 5 (36)] will prove most promising. Control of the interaction strength via the use of the Feshbach resonance may also be considered (33).
Acknowledgments
This work was supported by U.S. Department of Energy Grant FG05-86ER45234 and National Science Foundation Grant DMR-0205328.
Footnotes
- †To whom correspondence may be addressed. E-mail: constantine.yannouleas{at}physics.gatech.edu or uzi.landman{at}physics.gatech.edu
-
Author contributions: C.Y. and U.L. designed research, performed research, and wrote the paper.
-
Conflict of interest statement: No conflicts declared.
-
This paper was submitted directly (Track II) to the PNAS office.
- Abbreviations:
- BEC,
- Bose–Einstein condensate;
- EXD,
- exact diagonalization;
- GHL,
- generalized Heitler–London;
- GP,
- Gross–Pitaevskii;
- QD,
- quantum dot;
- ST,
- singlet–triplet;
- UHF,
- unrestricted Hartree–Fock;
- WM,
- Wigner molecule.
Abbreviations:
- © 2006 by The National Academy of Sciences of the USA