text:
# Local absorption spectra of artificial atoms and molecules

## I Introduction

In recent years much attention is being devoted to the properties of semiconductor quantum dots (QDs). In these systems, carriers are subject to a confining potential in all spatial directions, giving rise to a discrete energy spectrum (“artificial atoms”) and to novel phenomena of interest for fundamental physics as well as for applications to electronic and optoelectronic devices . The extension and the shape of the QD confining potential varies, depending on the nanostructure fabrication technique: The dots that are studied most extensively by optical methods are induced by quantum well (QW) thickness fluctuations , or obtained by spontaneous island formation in strained layer epitaxy , self-organized growth on patterned substrates , stressor-induced QW potential modulation , cleaved edge overgrowth , as well as chemical self-aggregation techniques . The resulting confinement lengths fall in a wide range between 1$`\mu `$m and 10 nm.

In spite of the continuing progress, all the available fabrication approaches still suffer from the effects of inhomogeneity and dispersion in the dot size, which lead to large linewidths when optical experiments are performed on large QD ensembles. A major advancement in the field has come from different types of local optical experiments, that allow the investigation of individual quantum dots thus avoiding inhomogeneous broadening .

Among local spectroscopies, the approaches based on scanning near-field optical microscopy (SNOM)  are especially interesting as they bring the spatial resolution well below the diffraction limit of light: With the development of small-aperture optical fiber probes, sub-wavelength resolutions were achieved ($`\lambda /8\lambda /5`$  or $`\lambda /40`$ ) and the first applications to nanostructures became possible . As the resolution increases, local optical techniques in principle allow direct access to the space and energy distribution of quantum states within the dot. This opens, however, a number of questions regarding the interpretation of these experiments, that were often neglected in the past.

First of all, for spatially inhomogeneous electromagnetic (EM) fields it is no longer possible to define and measure an absorption coefficient that locally relates the absorbed power density with the light intensity (since the susceptibility $`\chi (𝐫,𝐫^{})`$ cannot be approximated by a local tensor). In the linear regime, a local absorption coefficient can still be defined, which is however a complicated function that depends on the specific EM field distribution . The interpretation of near-field spectra therefore requires calculations based on a reasonable assumption for the profile of the EM field.

Secondly, the quantum states that are actually probed are few-particle states of the interacting electrons and holes photoexcited in the dot. Even in the linear regime, excitonic effects are known to dominate the optical spectra of dots since Coulomb interactions are strongly enhanced by the three-dimensional confinement. Near-field spectra probe exciton wavefunctions, and their spatial coherence and overlap with the em-field profile will determine the local absorption .

In this paper, we show how the above phenomena affect local spectra of QDs, with special attention to the case of coupled dots (“artificial molecules”) where carriers interact across the barrier via tunneling and/or Coulomb matrix elements . Indeed, the optical properties of coupled dots are currently of great interest not only in view of the unavoidable inter-dot interactions occurring in real samples with dense QD packing, but also in view of their relevance for designing novel devices including gates for possible solid-state implementations of quantum information processing .

We will show that the relative phase of the exciton wavefunction in adjacent coupled dots (or in different regions of the same dot) can induce dramatic changes in the selection rules with respect to far-field spectra: A realistic prediction of these effects require accurate calculations taking into account quantum confinement as well as Coulomb interactions. Our theoretical scheme is especially designed to allow a realistic description of the quantum states of the interacting electron and holes photoexcited in the linear regime. In this respect we improve drastically over previous approaches, which generally focused on a more detailed treatment of the EM field distributions .

Our theoretical framework for dots is summarized in Section II, while Sections III and IV discuss our results and conclusions for single and coupled dots.

## II Theory

In this section, we summarize our theoretical approach for computing local absorption spectra for semiconductor QDs. We first show in Sec. II A how to compute the single-particle eigenstates for electrons and holes subjected to a three-dimensional confinement potential. These single-particle states are then used in Sec. II B for the calculation of electron-hole (i.e., optical) excitations. In analogy to semiconductor systems of higher dimensionality, we shall refer to these excitations as excitons; the properties of such excitons, however, are not only governed by the attractive electron-hole Coulomb interaction but in addition by the strong quantum confinement. Finally, we use in Sec. II C the above ingredients to derive the equations needed for the calculation of local optical absorption spectra.

### A Single-particle states

In semiconductor QDs, carriers are confined in all three space directions. To simplify our analysis, we assume that a suitable parameterization of the dot confinement potential is known (e.g., from experiment) and that the confinement potential varies sufficiently slow on the length scale of the lattice constant. We thus shall make use of the envelope-function approach ; moreover, since the energy region of our present concern is relatively close to the semiconductor band gap, we describe the material band structure in terms of a single electron and hole band within the usual effective-mass approximation. More specifically, the envelope-function equation for single electrons and holes reads:

$$\left(\frac{\mathrm{}^2^2}{2m_{e,h}}+V_c^{e,h}(𝐫)\right)\varphi _\mu ^{e,h}(𝐫)=ϵ_\mu ^{e,h}\varphi _\mu ^{e,h}(𝐫),$$
(1)
where $`m_e`$ ($`m_h`$) is the effective mass and $`V_c^e`$ ($`V_c^h`$) is the confinement potential energy for electrons (holes). Following our approach developed earlier , we numerically solve Eq. (1) for arbitrary confinement potentials by use of a plane-wave expansion with periodic boundary conditions (see Appendix A).

### B Exciton states

When the dot structure is perturbed by an external light field (e.g., laser), electron-hole pairs are created which propagate in the presence of the mutual Coulomb interaction and of the dot confinement potential. Within the present paper, we shall restrict ourselves to the linear optical response, i.e., the dynamics of a single electron-hole pair. Then, the exciton dynamics is described by the electron-hole wavefunction $`\mathrm{\Psi }(𝐫_e,𝐫_h)`$, with the squared modulus being the probability of finding the electron at position $`𝐫_e`$ when the hole is at position $`𝐫_h`$.

If we expand the electron-hole (“exciton”) eigenfunction in terms of single-particle states viz.:

$$\mathrm{\Psi }^\lambda (𝐫_e,𝐫_h)=\underset{\mu \nu }{}\varphi _\mu ^e(𝐫_e)\mathrm{\Psi }_{\mu \nu }^\lambda \varphi _\nu ^h(𝐫_h),$$
(2)
we obtain the excitonic eigenvalue problem :

$$(ϵ_\mu ^e+ϵ_\nu ^h)\mathrm{\Psi }_{\mu \nu }^\lambda +\underset{\mu ^{}\nu ^{}}{}V_{\mu \mu ^{},\nu \nu ^{}}^{eh}\mathrm{\Psi }_{\mu ^{}\nu ^{}}^\lambda =E_\lambda \mathrm{\Psi }_{\mu \nu }^\lambda .$$
(3)
As will be shown in the following, the exciton spectrum $`E_\lambda `$ directly provides the optical transition energies whereas the excitonic wavefunctions $`\mathrm{\Psi }^\lambda `$ determine the oscillator strengths of the corresponding transitions. In Eq. (3) we have introduced the electron-hole Coulomb matrix elements :

$$V_{\mu \mu ^{},\nu \nu ^{}}^{eh}=e^2𝑑𝐫_e𝑑𝐫_h\frac{\varphi _{\mu }^{e}{}_{}{}^{}(𝐫_e)\varphi _\mu ^{}^e(𝐫_e)\varphi _{\nu }^{h}{}_{}{}^{}(𝐫_h)\varphi _\nu ^{}^h(𝐫_h)}{\kappa _o|𝐫_e𝐫_h|},$$
(4)
where $`e`$ is the elementary charge and $`\kappa _o`$ is the static dielectric constant of the bulk semiconductor (note that in Eq. (4) we have not considered the electron-hole exchange interaction). Within our computational approach, we consider in Eq. (3) typically a basis of 12 states for electrons and holes, respectively, and obtain the excitonic eigenfunctions by direct diagonalization of the Hamiltonian matrix.

### C Local optical absorption

When the semiconductor nanostructure is excited by a local near-field probe, the total absorbed power $`\alpha (\omega )`$ at a given frequency $`\omega `$ is proportional to $`𝑑𝐫_\omega (𝐫)P(𝐫,\omega )`$, where $`_\omega (𝐫)`$ is the electro-magnetic field distribution of the near-field probe. Within linear response, the induced interband polarization $`P(𝐫,\omega )`$ is related to $`_\omega (𝐫)`$ through:

$$P(𝐫,\omega )=𝑑𝐫^{}\chi (𝐫,𝐫^{};\omega )_\omega (𝐫^{}),$$
(5)
where the non-local electrical susceptibility $`\chi (𝐫,𝐫^{};\omega )`$ can be expressed in terms of the excitonic eigenenergies and eigenfunctions :

$$\chi (𝐫,𝐫^{};\omega )=\mu _o^2\underset{\lambda }{}\frac{\mathrm{\Psi }^\lambda (𝐫,𝐫)\mathrm{\Psi }_{}^{\lambda }{}_{}{}^{}(𝐫^{},𝐫^{})}{E_\lambda \mathrm{}\omega i\gamma }.$$
(6)
Here, $`\mu _o`$ is the dipole-matrix element of the bulk semiconductor, while we have introduced a small damping constant $`\gamma `$ accounting for the finite lifetime of exciton states due to environment coupling (e.g., phonons). To derive our final expression, it turns out to be convenient to consider for the elctromagnetic field distribution a given profile $`\xi `$ centered around the beam position $`𝐑`$, i.e., $`_\omega (𝐫)=_\omega \xi (𝐫𝐑)`$. Then, the local spectrum for a given tip position $`𝐑`$ can be expressed in the form  (see also Appendix A):

$$\alpha _\xi (𝐑,\omega )\mathrm{}\underset{\lambda }{}\frac{\alpha _\xi ^\lambda (𝐑)}{E_\lambda i\gamma \mathrm{}\omega },$$
(7)
where

$$\alpha _\xi ^\lambda (𝐑)=\left|𝑑𝐫\mathrm{\Psi }^\lambda (𝐫,𝐫)\xi (𝐫𝐑)\right|^2.$$
(8)
Two limiting cases can be identified. For a spatially homogeneous electromagnetic field (far-field), the oscillator strength $`\alpha _\xi ^\lambda `$ is given by the spatial average of the excitonic wavefunction, i.e., $`\alpha _\xi ^\lambda (𝐑)=|𝑑𝐫\mathrm{\Psi }^\lambda (𝐫,𝐫)|^2`$. In the opposite (and hypothetical) limit of an infinitely narrow probe, $`\xi (𝐫𝐑)=\delta (𝐫𝐑)`$, one is probing the local value of the exciton wavefunction, i.e., $`\alpha _\xi ^\lambda (𝐑)=|\mathrm{\Psi }^\lambda (𝐑,𝐑)|^2`$. Finally, within the intermediate regime of a narrow but finite probe, $`\mathrm{\Psi }^\lambda (𝐫,𝐫)`$ is averaged over a region which is determined by the spatial extension of the light beam; therefore, excitonic transitions which are optically forbidden in the far-field may become visible in the near-field.

## III Results

In the following sections we consider the interaction of the EM field with excitonic states of single and double QDs; for the latter system, we focus particularly on the transition between two isolated QDs and a the “artificial molecule”, where the electronic states of two QDs are strongly overlapping.

### A Single-particle states

We shall consider a prototypical QD confinement which is composed of a 2D harmonic potential in the $`(x,y)`$-plane and a rectangular quantum well along $`z`$; such confinement potentials have been demonstrated to be a good approximation for self-assembled quantum dots formed by strained-layer epitaxy. We focus on cases where the $`z`$-confinement is stronger than the $`(x,y)`$ one, so that the confinement potential can be written as

$$V_c^{e,h}(x,y,z)=V_{}^{e,h}(x,y)+V_o^{e,h}\theta (|z|\frac{z_o}{2}),$$
(9)
where $`z_o`$ is the width of the quantum-well and $`V_o^{e,h}`$ the band offsets for electrons and holes, respectively. For a single dot the in-plane confinement potential $`V_{}^{e,h}(x,y)`$ is of the form:

$$V_{}^{e,h}(x,y)=\frac{1}{2}𝒦_{e,h}(x^2+y^2),$$
(10)
while for two dots (i.e., double dot) separated by the distance $`d`$:

$$V_{}^{e,h}(x,y)=\{\begin{array}{cc}\frac{1}{2}𝒦_{e,h}\left((|x|\frac{d}{2})^2+y^2\right)\hfill & \text{for }|x|>\frac{d}{4}\hfill \\ \frac{1}{2}𝒦_{e,h}\left((\frac{d^2}{8}x^2)+y^2\right)\hfill & \text{otherwise}\hfill \end{array}$$
(11)
with $`𝒦_{e,h}=m_{e,h}(\omega _o^{e,h})^2`$, and $`\mathrm{}\omega _o^{e,h}`$ the level splittings of the in-plane harmonic potential. The shape of the double-dot potential has been obtained by matching the parabolas with opposite curvature, such that the potential is continuous and smooth at $`x=\pm \frac{d}{4}`$; the shape of the resulting potential along the $`x`$\- direction is shown, for selected inter-dot distances $`d`$, in Figs. 1(b), 1(d). Material and dot parameters which are used in this paper are listed in Table I; with this choice of parameters, electron and hole wavefunctions have approximately the same lateral extension, and the QW-induced intersubband splittings are much larger than $`\mathrm{}\omega _o^e`$ and $`\mathrm{}\omega _o^h`$.

With our choice of the confinement potential, Eq. (refeq:total.confinement), the single-particle energies of a QD are $`E_{\text{QD}}=E_{\text{QW}}+E_{\text{harm}}`$, where $`E_{\text{QW}}`$ is the confinement energy of the QW along $`z`$ and $`E_{\text{harm}}`$ is the confienment energy of the 2D single- or double-harmonic potential. Single-particle energies and envelope functions have been computed numerically within a plane-wave scheme. However, for a single QD the 2D eigenstates can be found analytically and are the well-known “Fock-Darwin” states  (we stress, however, that the extension in the $`z`$-direction is of crucial importance for the calculation of the Coulomb matrix elements and the optical properties, and unavoidably has to be taken into account in any realistic calculation; see also discussion in Ref. ). For such states $`E_{\text{harm}}=(n+1)\mathrm{}\omega _o^{e,h}`$, where $`n=0,1,\mathrm{}`$ is the principal quantum number, and each level is $`(n+1)`$-fold degenerate; in Table LABEL:table:fock.darwin we summarize for convenience some properties of these “Fock-Darwin” states .

Figure 1 shows the calculated single-particle energies for electrons and holes for the more complex case of a double QD with the confinement potential given in Eq. (11) and with parameters listed in Table I. The lower panels (Figs. 1(b,d)) show the confinement potentials for electrons and holes at selected inter-dot distances. Obviously, for large dot separations $`d`$ ($`d60`$ nm) the system can be well approximated by two separate QDs; in this regime the equidistance of the excited states and the correct degeneracy of the Fock-Darwin states is obtained. When $`d`$ is small enough that carriers have sufficient energy to overcome (or tunnel through) the barrier between the two dots, the degeneracy is removed, and the energy levels have a non-monotonous behaviour which reflects the transition from two separated carrier systems to a single one, and is similar to the one found, e.g. for coupled QWs . For the smallest dot distances the double-dot potential merges into a single-dot potential, and the Fock-Darwin states of a single dot are recovered.

### B Role of the Coulomb correlation in the far-field spectra

Before turning to the analysis of near-field spectra, we shortly discuss the limiting case of very broad EM field distribution (far-field spectra). This discussion allows us to elucidate the role of the electron-hole Coulomb correlation, particularly in the transition from two separate “artificial atoms” to an “artificial molecule”.

Far-field spectra can be obtained in the formalism of Sec. II C with a spatially homogeneous electromagnetic field distribution probe $`\xi `$):

$$\alpha _{\xi (𝐫)=\mathrm{const}}^\lambda =\left|𝑑𝐫\mathrm{\Psi }^\lambda (𝐫,𝐫)\right|^2.$$
(12)
Figs. 2 and 3 show the calculated far-field spectra for a double QD as a function of the dot distance $`d`$. We first concentrate on the calculations where Coulomb correlations were artificially set to zero (Figs. 2(a) and 3(a)): Because of symmetry, only a small fraction of all possible electron-hole transitions is visible; from Eq. (12) and using $`𝑑\phi \mathrm{exp}i(m_e+m_h)\phi \delta _{m_e,m_h}`$ we obtain that optical transitions are only allowed between electron and hole single-particle states with opposite angular momentum. Indeed, for large distances (uncoupled QDs) only three strong absorption peaks are observed, with an energy splitting of approximately $`\mathrm{}\omega _o^e+\mathrm{}\omega _o^h`$; the intensity of the peaks increases with energy (with ratio $`1:2:3`$). These can be attributed to transitions between electrons and holes single-particle states (see Table LABEL:table:fock.darwin) of the $`1s`$ symmetry (peak at $`70`$ meV), the $`1p`$ symmetry (peak at $`95`$ meV), and the $`1d`$ and $`2s`$ symmetries (peak at $`120`$ meV).

When symmetry is reduced, either because of an asymmetric confinement potential or by the presence of an external inhomogeneous EM-field (as will be discussed later), the selection rules noted above are relaxed. Indeed, when $`d`$ is reduced and the two QDs begin to interact, the calculated spectra show a much richer structure, as shown in Fig. 2(a) and 3(a), reflecting the reduction of built-in symmetry. Obviously, when $`d0`$, the usual selection rules of a single-QD are recovered.

When Coulomb interaction is included, inspection of the exciton wavefunctions $`\mathrm{\Psi }_{\mu \nu }^\lambda `$ (obtained from the solutions of Eq. (3)) shows that a number of different single-particle transitions contributes to each excitonic state . Coulomb interaction affects the optical spectra \[see Figs. 3(a) and 3(b)\] in several ways. Firstly, because of the attractive electron-hole interaction leading to the “bound” excitonic states, we observe a red-shift of the peaks; peak separation is barely affected, however, at least for large $`d`$. Secondly, we observe a redistribution of the oscillator strength; for isolated QDs the three main peaks are of similar height. In general, oscillator strength is transferred from higher to lower peaks. This effect is particularly strong, e.g., in the doublet which splits from the lowest peak when the two QDs approach; contrary to the uncorrelated case, the heighest partner is very weak. Finally, Coulomb interaction is responsible for the appearance of additional lines (see, e.g., for $`d=70`$ nm the peak at 70 meV). While the first two effects (red-shift and transfer of oscillator strength) are similar to what is found in the absorption spectra of semiconductor quantum wires, and thus can be considered as a general fingerprint of Coulomb correlations in the optical properties of semiconductor nanostructures, the origin of the additional peaks is best discussed in connection with the calculated near-field optical spectra and, therefore, is postponed to the next section.

### C Optical near-field spectra

In this section we discuss the local absorption spectra of single and coupled QDs. Because of the narrow well width of the dot confinement potential (see Table I), the EM profile of the near-field probe along $`z`$ has an only minor influence on the results, and we use:

$$\xi (x,y,z)\mathrm{exp}\left(\frac{x^2+y^2}{2\sigma ^2}\right).$$
(13)
The spatial resolution of the electromagnetic field distribution of Eq. (13) is then approximately given by the full-width at half maximum (FWHM) of the Gaussian (i.e., $`2\sqrt{2\mathrm{ln}2}\sigma 2.35\sigma `$). Since the Gaussian acts as an envelope on $`\mathrm{\Psi }^\lambda `$, in the intermediate regime of a narrow but finite $`\sigma `$ the spatial average only extends over the region where the Gaussian is non-vanishing.

Since the extension of the quantum states under investigation is of the order of a few tens of nano-meters (see also Figs. 5 and 8, to be discussed below), in our calculations we consider three different regimes of spatial resolution: i) a regime where the FWHM is much larger than the extension of the quantum states (as a characteristic value we use $`\sigma =50`$ nm); ii) a regime where the FWHM is comparable to the extension of the relevant quantum states (we use $`\sigma =10`$ nm); iii) a regime with an extremely narrow probe beam (we use $`\sigma =0.1`$ nm). Calculations performed in this latter (unphysical) regime are used for illustrative purposes to obtain a “cartography” of the exciton wavefunction, as discussed at the end of Sec. II C. We finally notice that the excitonic Bohr radius $`12`$ nm for GaAs.

#### 1 Single quantum dot

In Fig. 4 we report the calculated local absorption spectra $`\alpha _\xi (X,\mathrm{}\omega )`$ for a single QD as a function of the tip position. The tip is swept along one direction, passing through the center of the QD.

In Figs. 4(a–c) we show the calculated spectra neglecting Coulomb interaction. For the highest spatial resolution \[Fig. 4(a)\], the local absorption at photon energy $`E_\lambda `$ is proportional to $`𝑑z\mathrm{\Psi }^\lambda (𝐫,𝐫)|_{y=0}`$. Given the energy splitting $`\mathrm{}\omega _o^h=3.5`$ meV for holes and $`\mathrm{}\omega _o^e=20`$ meV for electrons, we can attribute the triplet of peaks at $`70`$ meV to the single-particle transitions involving the $`1s`$ state of electrons and the $`1s`$, $`1p`$, and $`(2s,1d)`$ states (in order of increasing energy) of holes (see also Table LABEL:table:fock.darwin); analogously, the triplet at $`90`$ meV is attributed to the transitions involving the $`1p`$ state of electrons and the $`1s`$, $`1p`$, and $`(2s,1d)`$ hole states; indeed, in Fig. 4(a) the localization of the absorption peaks is suggestive of the $`s`$-, $`p`$\- or $`d`$-type symmetry of the corresponding Fock-Darwin states. These features are still present at the intermediate resolution \[Fig. 4(b)\], but disappear at the opposite limit of a broad probe \[Figs. 4(c)\]. This is expected, since, when a localized EM-field is present, the symmetry of the whole system (nanostructure+EM-field) is lower than that of the nanostructure (except when the probe is centered in the symmetry center of the structure), and far-field selection rules are relaxed. When the probe is broadened, however, the built-in symmetry of the structure is recovered, and optical far-field selection rules (i.e., optical transitions only between electron and hole states with opposite angular momentum $`m`$) apply; therefore, the spectra are almost identical to those of two separated dots in far-field spectroscopy, already discussed in Figs. 2 and 3.

When we compare Figs. 4(a–c) with Figs. 4(d–f), we find that Coulomb interaction induces several effects which are expected on the basis of the discussion of the far-field spectra. In particular, we find i) an almost rigid redshift of the spectra; ii) a transfer of oscillator strength from transitions at higher energies to those at lower energies; iii) the appearance of new features in the optical spectra. To discuss the origin of these new optical features caused by Coulomb interactions, let us consider, e.g., the optical peaks at photon energy $`65`$ meV (Figs. 4(d–f)): They are quite strong at $`\sigma =0.1`$ nm (Fig. 4(d)), almost disappear at $`\sigma =10`$ nm (Fig. 4(e)), and are visible again in the far field limit (Fig. 4(f)). Such a behaviour is rather unexpected and noticeably differs from that of other transitions, which —with increasing $`\sigma `$— either remain strong or gradually disappear due to symmetry reasons, as discussed above. To investigate the origin of this non-monotonic dependence, in the following we analyze the three excitons within the corresponding energy range. Figure 5 shows a contour plot of the respective exciton wavefunction $`\mathrm{\Psi }^\lambda (𝐫,𝐫)|_{z=0}`$. Apparently, in Fig. 5(a) the exciton has $`s`$-type symmetry, whereas the other two electron-hole states have $`p`$-type symmetry. (Because of the periodicity box used in our calculations, the two-fold degenerate $`p`$-type exciton wavefunctions have cartesian rather than cylinder symmetry; note that, since the presence of the near-field tip destroys the cylinder symmetry, the wavefunctions shown in Fig. 5 indeed form a natural basis; see also Table LABEL:table:fock.darwin). Next, we note that the average $`𝑑𝐫\mathrm{\Psi }^\lambda (𝐫,𝐫)`$ of the $`p`$-type exciton wavefunctions is zero. Since with increasing $`\sigma `$ the radius within which the exciton eigenfunctions $`\mathrm{\Psi }^\lambda `$ are averaged increases, we expect for these $`p`$-type functions with increasing $`\sigma `$ a monotonically decreasing behavior. The exciton shown in Fig. 5(a), on the other hand, has a non-zero average and is therefore visible in both the optical far- and near-field. A closer inspection of the exciton wavefunction $`\mathrm{\Psi }_{\mu \nu }^\lambda `$ reveals that the largest contribution stems from the transition between the $`1s`$ state of electrons and the $`2s`$ state of holes, but there is also a noticeable contribution from the $`1s`$-$`1s`$ and $`1p`$-$`1p`$ electron-hole transitions. Indeed, only the latter contributions couple in the far-field to the light field. In the regime of finite resolution, there is an optimal cancellation when the FWHM of the EM near field (i.e. $`2\sqrt{2\mathrm{ln}2}\sigma 2.35\sigma `$) becomes equal to the Bohr radius. This is clearly depicted in Fig. 6, where, in order to facilitate our discussion, we have introduced the quantity $`I_\xi ^\lambda 𝑑𝐑\alpha _\xi ^\lambda (𝐑)`$, which provides a measure of the relative contribution of each exciton to the absorption spectra. Fig. 6 shows $`I_\xi ^\lambda `$ for the three excitons (shown Fig. 5) within the energy region of 65 meV: We observe that with increasing $`\sigma `$, the $`p`$-type functions (open circles) indeed vanish monotonically, whereas for the $`s`$-type exciton (full circles) there exists an optimal cancellation when the FWHM of the EM field distribution becomes approximately equal to the Bohr radius. In spite of the specific carrier states of a single parabolic QD, we expect that such non-monotonic behavior appears quite generally in semiconductor nanostructures where carrier states are confined on a length scale comparable to the Bohr radius, and thus provides a striking fingerprint of Coulomb correlations in the optical near-field spectra (we find similar behavior in our calculations for the near-field spectra of coupled QDs discussed below).

#### 2 Double quantum dot

In Figure 7 we show the calculated local absorption spectra $`\alpha _\xi (X,\mathrm{}\omega )`$ for a double QD for selected values of the interdot distance and $`\sigma `$. The tip of the probe is swept along the direction which passes through the centers of the two QDs.

Let us first concentrate on the results with $`\sigma =0.1`$ nm and with the Coulomb interaction taken into account (Figs. 7(d,g,j)). With decreasing interdot distance we observe the transition from a system where the energetically lowest exciton states are almost localized in the spatially separated minima of the two dots, to a system where the electron-hole states extend over the whole nanostructure. Here, the $`s`$-like ground-state excitons of Fig. 7(j) split up into a “bonding” and an “anti-bonding” state (Fig. 7(d)). By comparing Figs. 7(d) and 7(f), we find that in the optical far-field only the symmetric ground state exciton couples to the light field.

Next, we discuss the optical features at the photon energy of $`70`$ meV for $`d=40`$ nm. As in the case of the single dot, these features show a non-monotonic dependence on the probe width. As can be inferred from the calculations with $`\sigma =0.1`$ nm, there are several excitonic states contributing to the spectral features in this energy range; Fig. 8 shows the excitonic wavefunction of two states out of the six states with $`E_\lambda 70`$ meV for for $`d=40`$ nm; it can be inferred that for a spatial resolution of the near-field probe comparable to the excitonic Bohr radius ($`12`$ nm) there is again an optimal cancellation. This is a remarkable finding, because it clearly demonstrates that such a behaviour indeed is a general characteristics of semiconductor nanostructures, and does not depend on peculiar symmetries of the confining potential.

## IV Summary and Conclusions

We have analyzed theoretically the interaction between a model near-field probe and a zero-dimensional heterostructure: Quantum confinement of the electron and hole states, as well as their Coulomb interaction in the linear regime are fully included in our description.

We have specifically considered single and coupled semiconductor quantum dots, and shown that absorption is strongly influenced by the spatial interference of the exciton wavefunctions, which depends on the spatial extension of the light beam. As a consequence, near-field experiments on quantum dots are predicted to display unexpected spectral features whose dependence on spatial resolution is highly non-trivial.

When combined with an appropriate choice of the EM field distribution, our approach provides the necessary tool for interpretation of near-field absorption spectra of quantum dots as the spatial resolution of experiments becomes comparable with the Bohr radius of the exciton in the nanostructure.

###### Acknowledgements.

We thank Fausto Rossi for most stimulating discussions. This work was supported in part by INFM through PRA-99-SSQI, and by the EC under the TMR Network “Ultrafast Quantum Optoelectronics” and the IST programme “SQID”. U.H. acknowledges support by the EC through a TMR Marie Curie Grant.

## A Plane-wave approach

In this Appendix we discuss details of our numerical solution schemes based on a plane-wave expansion. Following our approach developed earlier , we consider the problem of a single or double QD which is located inside a box with periodic boundary conditions, where the boxsize is chosen sufficiently large to avoid interactions with “neighbour” dots. As a complete set of functions, inside the periodicity box we use a plane-wave basis, $`|𝐤`$, with:

$$k_\alpha =\frac{2\pi n_\alpha }{L_\alpha },n_\alpha 𝒵,\alpha =x,y,z.$$
(A1)
Here $`L_\alpha `$ denotes the sizes of the periodicity box (we use the same box for electrons and holes). We next expand the single-particle wavefunctions for electrons and holes within the plane-wave basis:

$$\stackrel{~}{\varphi }_{\mu ,𝐤}^{e,h}=\mathrm{\Omega }^1𝑑𝐫e^{i𝐤𝐫}\varphi _\mu ^{e,h}(𝐫),$$
(A2)
with $`\mathrm{\Omega }`$ the volume of the periodicity box. The envelope-function equation (1) is then transformed to:

$$\underset{𝐤^{}}{}\left(\frac{\mathrm{}^2𝐤^2}{2m_{e,h}}\delta _{\mathrm{𝐤𝐤}^{}}+\stackrel{~}{V}_{c,𝐤𝐤^{}}^{e,h}\right)\stackrel{~}{\varphi }_{\mu ,𝐤^{}}^{e,h}=ϵ_\mu ^{e,h}\stackrel{~}{\varphi }_{\mu ,𝐤}^{e,h},$$
(A3)
which can be solved by standard diagonalization techniques. To keep the numerics tractable, only wavevectors smaller than a given cut-off wavevector are considered (typically 2000–3000 wavevectors). In our computational approach, we perform the Fourier transform of the confinement potential by storing $`V_c^{e,h}(𝐫)`$ on an appropriate grid (with a typical number of 30 points along each direction), and approximating within each cube $`V_c^{e,h}(𝐫)`$ by its average value.

In the calculation of the near-field spectra, we define the electron-hole index $`l=(\mu ,\nu )`$. Then:

$$\mathrm{\Psi }^\lambda (𝐫,𝐫)=\underset{l}{}\mathrm{\Psi }_l^\lambda \varphi _{\mu _l}^e(𝐫)\varphi _{\nu _l}^h(𝐫),$$
(A4)
and we obtain for $`\alpha _\xi ^\lambda (𝐑)`$ of Eq. (8) the final result:

$$\alpha _\xi ^\lambda (𝐑)=|\underset{l}{}\mathrm{\Psi }_l^\lambda \underset{𝐤,𝐤^{}}{}\stackrel{~}{\xi }_{𝐤+𝐤^{}}(𝐑)\stackrel{~}{\varphi }_{\mu _l,𝐤}^e\stackrel{~}{\varphi }_{\nu _l,𝐤^{}}^h|^2,$$
(A5)
with $`\stackrel{~}{\xi }_𝐤(𝐑)=\mathrm{\Omega }^1𝑑𝐫\xi (𝐫)e^{i𝐤(𝐫+𝐑)}`$.


format: text

text:
# Many-Body Theory of Dilute Bose-Einstein Condensates with Internal Degrees of Freedom

## Abstract

The Bogoliubov theory of weakly interacting bosons is generalized to Bose-Einstein condensates with internal degrees of freedom so that a single effective Hamiltonian produces various many-body ground states or metastable spin domains and the corresponding collective modes on an equal footing.

The internal degrees of freedom of alkali Bose-Einstein condensates (BECs) arise primarily from electronic spin and are hence by far more amenable to manipulation than those of superfluid 3He. The realization of “spinor” BECs  has motivated many researchers to study the magnetism of superfluid vapors . The excitation spectra of spinor BECs have been examined by Ohmi and Machida  and by Ho  using mean-field theory. The many-body ground states were investigated by Law et al.  and by Koashi and Ueda  and Ho and Yip  in the absence and presence, respectively, of a magnetic field; in the latter case, it is predicted that the magnetic sublevel $`m=0`$ of a spin-1 antiferromagnetic BEC becomes populated as the magnetic field decreases. The main aim of this Letter is to present a unified theory that derives the various many-body states of BEC and the corresponding excitation spectra on an equal footing. Previous work  has been unable to address this problem, as it has been assumed that the internal degrees of freedom are independent of other degrees of freedom such as the density and spin waves.

We consider a system of $`N`$ spin-1 identical bosons interacting via a binary, contact interaction, subject to a uniform magnetic field $`𝐁`$. The Hamiltonian of the system is given by

$`\widehat{H}`$ $`=`$ $`{\displaystyle }d𝐫[\widehat{\mathrm{\Psi }}_\alpha ^{}({\displaystyle \frac{\mathrm{}^2}{2M}}^2+U(𝐫))\widehat{\mathrm{\Psi }}_\alpha p\widehat{\mathrm{\Psi }}_\alpha ^{}f_{\alpha \beta }^z\widehat{\mathrm{\Psi }}_\beta `$ (2)
$`+{\displaystyle \frac{c_0}{2}}\widehat{\mathrm{\Psi }}_\alpha ^{}\widehat{\mathrm{\Psi }}_\beta ^{}\widehat{\mathrm{\Psi }}_\beta \widehat{\mathrm{\Psi }}_\alpha +{\displaystyle \frac{c_1}{2}}\widehat{\mathrm{\Psi }}_\alpha ^{}\widehat{\mathrm{\Psi }}_\alpha ^{}^{}𝐟_{\alpha \beta }𝐟_{\alpha ^{}\beta ^{}}\widehat{\mathrm{\Psi }}_\beta ^{}\widehat{\mathrm{\Psi }}_\beta ],`$
where $`p(>0)`$ is the product of the gyromagnetic ratio and the external magnetic field $`B`$ applied in the z-direction; $`𝐟_{\alpha \beta }=(f_{\alpha \beta }^x,f_{\alpha \beta }^y,f_{\alpha \beta }^z)`$ represents the components ($`\alpha ,\beta =1,0,1`$) of spin-1 matrices; and $`\widehat{\mathrm{\Psi }}_\alpha (𝐫)`$ denotes the field operator that annihilates a boson in magnetic sublevel $`\alpha `$ at position $`𝐫`$. The quadratic Zeeman effect will be taken into account when it becomes relevant. The coupling constants $`c_0`$ and $`c_1`$ are related to the spin-singlet and triplet s-wave scattering lengths $`a_s`$ and $`a_t`$ by $`c_0=4\pi \mathrm{}^2(a_s+2a_t)/3M`$ and $`c_1=4\pi \mathrm{}^2(a_ta_s)/3M`$, respectively. In Eq. (2) the repeated indicies are assumed to be summed.

When the system has translation invariance, it is convenient to expand the field operator as $`\widehat{\mathrm{\Psi }}_\alpha (𝐫)=(1/\sqrt{V})_𝐤\widehat{a}_{𝐤\alpha }e^{i\mathrm{𝐤𝐫}}`$, where $`V`$ is the volume of the system. Equation (2) then reduces to

$`\widehat{H}={\displaystyle \underset{𝐤\alpha }{}}(ϵ_𝐤p\alpha )\widehat{n}_{𝐤\alpha }+{\displaystyle \frac{1}{2V}}{\displaystyle \underset{𝐤}{}}:(c_0\widehat{\rho }_𝐤^{}\widehat{\rho }_𝐤+c_1\widehat{𝐒}_𝐤^{}\widehat{𝐒}_𝐤):,`$ (3)
where $`ϵ_𝐤\mathrm{}^2𝐤^2/2M`$, $`\widehat{n}_{𝐤\alpha }\widehat{a}_{𝐤\alpha }^{}\widehat{a}_{𝐤\alpha }`$, $`\widehat{\rho }_𝐤_𝐪\widehat{a}_{𝐪\alpha }^{}\widehat{a}_{𝐪+𝐤\alpha }`$, $`\widehat{𝐒}_𝐤_{𝐪\alpha \beta }𝐟_{\alpha \beta }\widehat{a}_{𝐪\alpha }^{}\widehat{a}_{𝐪+𝐤\beta }`$, and $`::`$ denotes normal ordering.

When the interparticle interactions are repulsive, BEC occurs in the $`𝐤=0`$ state. We therefore separate the $`𝐤=0`$ components $`\widehat{a}_{\mathrm{𝟎}\alpha }`$ in the Hamiltonian and keep the terms only up to second order in $`\widehat{a}_{𝐤\alpha }`$ ($`𝐤0`$). Because our interest is to derive various many-body states and the corresponding low-lying excitation spectra on an equal footing, we should not replace the operators $`\widehat{a}_{\mathrm{𝟎}\alpha }`$ by c-numbers. After some algebraic manipulation, we obtain

$`\widehat{H}^{\mathrm{eff}}={\displaystyle \frac{c_0+c_1}{2V}}N(N1){\displaystyle \frac{3c_1}{2V}}\widehat{A}^{}\widehat{A}p(\widehat{n}_{\mathrm{𝟎}1}\widehat{n}_{\mathrm{𝟎}1})`$ (4)
$`+{\displaystyle \underset{𝐤\mathrm{𝟎}}{}}{\displaystyle \underset{\alpha =0,\pm 1}{}}\left(ϵ_𝐤p\alpha {\displaystyle \frac{(1+|\alpha |)c_1}{V}}\widehat{n}_{\mathrm{𝟎}\alpha }\right)\widehat{n}_{𝐤\alpha }`$ (5)
$`+{\displaystyle \frac{1}{V}}{\displaystyle \underset{𝐤\mathrm{𝟎}}{}}\{{\displaystyle \frac{c_0}{2}}(\widehat{D}_𝐤\widehat{D}_𝐤+\widehat{D}_𝐤^{}\widehat{D}_𝐤)+{\displaystyle \frac{c_1}{2}}(\widehat{S}_𝐤\widehat{S}_𝐤+\widehat{S}_𝐤^{}\widehat{S}_𝐤)`$ (6)
$`+{\displaystyle \underset{\alpha =0,1}{}}c_1[\widehat{a}_{𝐤\alpha }^{}\widehat{a}_{𝐤\alpha }^{}\widehat{a}_{\mathrm{𝟎}1\alpha }\widehat{a}_{\mathrm{𝟎}\alpha 1}+\widehat{a}_{𝐤\alpha 1}^{}\widehat{a}_{𝐤\alpha }^{}\widehat{a}_{\mathrm{𝟎}\alpha 1}\widehat{a}_{\mathrm{𝟎}\alpha }`$ (7)
$`+(\widehat{a}_{𝐤\alpha }^{}\widehat{a}_{𝐤\alpha 1}+2\widehat{a}_{𝐤1\alpha }^{}\widehat{a}_{𝐤\alpha })\widehat{a}_{\mathrm{𝟎}\alpha 1}^{}\widehat{a}_{\mathrm{𝟎}\alpha }]+\mathrm{h}.\mathrm{c}.\}`$ (8)
where h.c. denotes the hermitian conjugates of all preceding terms in the curly brackets, $`\widehat{A}^{}(\widehat{a}_{\mathrm{𝟎}0}^22\widehat{a}_{\mathrm{𝟎}1}^{}\widehat{a}_{\mathrm{𝟎}1}^{})/\sqrt{3}`$ creates a pair of spin-singlet bosons when acting upon the vacuum, and $`\widehat{D}_𝐤_\alpha \widehat{a}_{0\alpha }^{}\widehat{a}_{𝐤\alpha }`$ and $`\widehat{S}_𝐤_\alpha \alpha \widehat{a}_{0\alpha }^{}\widehat{a}_{𝐤\alpha }`$ describe density-wave and spin-wave operators, respectively. The first three terms on the right-hand side of Eq. (8) were studied in Refs. ; the remaining terms show how the internal and spatial degrees of freedom couple to one another.

To set a reference frame for later discussions, let us recapitulate the main results of Refs.. The Hamiltonian discussed in these references reads essentially

$`\widehat{H}_0={\displaystyle \frac{3c_1}{2V}}\widehat{A}^{}\widehat{A}p(\widehat{n}_{\mathrm{𝟎}1}\widehat{n}_{\mathrm{𝟎}1}).`$ (9)
When $`c_1<0`$, the energy minimum is obtained by putting all particles in the state with $`𝐤=\mathrm{𝟎}`$ and $`\alpha =1`$:

$`|N_2=0,F=N,F_z=N={\displaystyle \frac{(\widehat{a}_{\mathrm{𝟎}1}^{})^N}{\sqrt{N!}}}|\mathrm{vac},`$ (10)
where $`F`$ and $`F_z`$ denote the total spin and its projection on the $`z`$ axis, and $`N_2`$ is given in terms of $`F`$ and the total number of atoms $`N`$ by $`N_2=(NF)/2`$. Because all spins are aligned in the same direction, the system is ferromagnetic. When $`c_1>0`$, the first term on the right-hand side of Eq. (9) energetically favors the spin-singlet (or antiferromagnetic) correlation, whereas the second one favors the parallel-spin configurations. For a given external magnetic field ($`p`$), the energy minimum is attained for the state

$`|N_2=(NF)/2,F,F_z=F_0`$ (11)
$``$ $`(\widehat{A}^{})^{N_2}(\widehat{F}_{})^{FF_z}(\widehat{a}_{\mathrm{𝟎}1}^{})^F|\mathrm{vac},`$ (12)
where $`F`$ denotes an integer that satisfies $`F1/2<p/c_1<F+3/2`$, and $`\widehat{F}_{}(f^xif^y)_{\alpha \beta }\widehat{a}_{\mathrm{𝟎}\alpha }^{}\widehat{a}_{\mathrm{𝟎}\beta }`$. In both Eqs. (10) and (12), the depletion of the condensate due to interactions is not taken into account.

When the components with $`𝐤=\mathrm{𝟎}`$ are macroscopically occupied, we may replace $`\widehat{a}_{\mathrm{𝟎}\alpha }`$ in Eq. (8) with the c-numbers $`\sqrt{N}\zeta _\alpha `$. To find the correct excitation spectra, it is crucial to take the depletion of the condensate into account. This is done by setting

$`{\displaystyle \underset{\alpha }{}}|\zeta _\alpha |^2=1{\displaystyle \frac{1}{N}}{\displaystyle \underset{𝐤\mathrm{𝟎},\alpha }{}}\widehat{n}_{𝐤\alpha },`$ (13)
where the last term describes the fraction of the depletion. The energy of the condensate is given by

$`E_0={\displaystyle \frac{c_1N^2}{2V}}|\zeta _0^22\zeta _1\zeta _1|^2pN(|\zeta _1|^2|\zeta _1|^2),`$ (14)
where $`\zeta _\alpha `$ are determined by requiring that $`E_0`$ be minimized subject to constraint (13).

For $`c_1<0`$, the minimum of $`E_0`$ is reached when

$`\zeta _0=\zeta _1=0,|\zeta _1|^2=1{\displaystyle \frac{1}{N}}{\displaystyle \underset{𝐤\mathrm{𝟎},\alpha }{}}\widehat{n}_{𝐤\alpha }.`$ (15)
The Hamiltonian (8) then reduces to

$`\widehat{H}^\mathrm{F}={\displaystyle \frac{c_0+c_1}{2V}}N(N1)pN`$ (16)
$`+{\displaystyle \underset{𝐤0}{}}[ϵ_𝐤\widehat{n}_{𝐤1}+{\displaystyle \frac{(c_0+c_1)N}{2V}}(\widehat{a}_{𝐤1}\widehat{a}_{𝐤1}+\widehat{n}_{𝐤1}+\mathrm{h}.\mathrm{c}.)`$ (17)
$`+(ϵ+p)\widehat{n}_{𝐤0}+(ϵ+2p{\displaystyle \frac{2c_1N}{V}})\widehat{n}_{𝐤1}].`$ (18)
Diagonalizing this Hamiltonian, we find the excitation spectra as

$`E_{𝐤,1}^\mathrm{F}`$ $`=`$ $`\sqrt{ϵ_\mathrm{k}[ϵ_\mathrm{k}+2(c_0+c_1)n]},`$ (19)
$`E_{𝐤,0}^\mathrm{F}`$ $`=`$ $`ϵ_𝐤+p,`$ (20)
$`E_{𝐤,1}^\mathrm{F}`$ $`=`$ $`ϵ_𝐤+2p2c_1n,`$ (21)
in agreement with Ref.. We may gain some insight into the nature of the many-body correlation by writing the many-body wave function in the coordinate representation:

$`\mathrm{\Psi }(𝐫_1\alpha _1,\mathrm{},𝐫_N\alpha _N)=\mathrm{vac}|\widehat{\mathrm{\Psi }}_{\alpha _1}(𝐫_1)\mathrm{}\widehat{\mathrm{\Psi }}_{\alpha _1}(𝐫_1)|\mathrm{\Phi }.`$ (22)
Here $`|\mathrm{\Phi }`$ is given in the Bogoliubov approximation by

$`|\mathrm{\Phi }\mathrm{exp}\left(\varphi _0\widehat{a}_{\mathrm{𝟎}1}^{}{\displaystyle \underset{k_x>0}{}}\nu _k\widehat{a}_{𝐤1}^{}\widehat{a}_{𝐤1}^{}\right)|\mathrm{vac},`$ (23)
where $`\varphi _0^2=N[1(8/3)\sqrt{na_\mathrm{t}^3/\pi }]`$ and $`\nu _k=1+c_k\sqrt{c_k(c_k+2)}`$ with $`c_kk^2/(8\pi a_\mathrm{t}n)`$. In the dilute limit, the many-body wave function becomes very small unless all $`\alpha _k`$’s are equal 1, and we find that

$`\mathrm{\Psi }(𝐫_11,\mathrm{},𝐫_N1){\displaystyle \frac{\mathrm{exp}\left(\frac{4}{9}(3\pi 8)N\sqrt{na_t^3/\pi }\right)}{(2\pi NV^{2N})^{1/4}}}`$ (24)
$`\times \mathrm{exp}\left({\displaystyle \underset{i<j}{}}{\displaystyle \frac{a_t}{r_{ij}}}e^{r_{ij}/\xi }\right),`$ (25)
where $`r_{ij}|𝐫_i𝐫_j|`$ and the terms that are of the order of $`1/N`$ are ignored. This result clearly shows that two bosons strongly repel each other when their distance becomes smaller than the healing length $`\xi (8\pi a_tn)^{1/2}`$.

For $`c_1>0`$, minimizing $`F=E_0\delta _\alpha |\zeta _\alpha |^2`$, where $`\delta `$ is a Lagrange multiplier, yields

$`|\zeta _0|\left(e^{i(\varphi _1+\varphi _12\varphi _0)}+{\displaystyle \frac{\delta }{\sqrt{\delta ^2\gamma ^2}}}\right)=0,`$ (26)
where $`\varphi _\alpha \mathrm{arg}\zeta _\alpha `$ and $`\gamma 2p/(c_1n)`$. In the presence of the external magnetic field (i.e., $`\gamma 0`$), the minimum of $`F`$ is attained when

$`\zeta _0=0,|\zeta _{\pm 1}|^2={\displaystyle \frac{1}{2}}\pm {\displaystyle \frac{\gamma }{4}}{\displaystyle \frac{1}{2N}}{\displaystyle \underset{𝐤0\alpha }{}}\widehat{n}_{𝐤,\alpha },`$ (27)
and the corresponding Hamiltonian reads

$`\widehat{H}_{\gamma 0}^{\mathrm{AF}}={\displaystyle \frac{c_0}{2V}}N(N1)+{\displaystyle \frac{c_1}{2V}}N\left({\displaystyle \frac{\gamma ^2N}{4}}1\right){\displaystyle \frac{p\gamma N}{2}}`$ (28)
$`{\displaystyle \underset{𝐤\mathrm{𝟎},\alpha }{}}(ϵ_𝐤+c_1n\delta _{\alpha ,0})\widehat{n}_{𝐤\alpha }+{\displaystyle \frac{1}{2V}}{\displaystyle \underset{𝐤0}{}}\{c_0(\widehat{D}_𝐤\widehat{D}_𝐤+\widehat{D}_𝐤^{}\widehat{D}_𝐤)`$ (29)
$`+c_1(\widehat{S}_𝐤\widehat{S}_𝐤+\widehat{S}_𝐤^{}\widehat{S}_𝐤+n|\zeta _1\zeta _1|\widehat{a}_{𝐤,0}\widehat{a}_{𝐤,0})+\mathrm{h}.\mathrm{c}.\}.`$ (30)
The $`\alpha =0`$ mode is decoupled and its dispersion relation is given by

$`E_{𝐤,0}^{\mathrm{AF}}=\sqrt{ϵ_\mathrm{k}^2+2c_1nϵ_\mathrm{k}+p^2}.`$ (31)
The $`\alpha =\pm 1`$ modes are coupled, and the dispersion relations of the coupled modes are given by

$`E_{𝐤,\pm }^{\mathrm{AF}}=\sqrt{ϵ_\mathrm{k}^2+2(c_0+c_1)nϵ_\mathrm{k}\pm ϵ_\mathrm{k}\sqrt{n^2(c_0c_1)^2+{\displaystyle \frac{4p^2c_0}{c_1}}}},`$ (32)
(33)
again in agreement with Ref. . Our theory thus reproduces all known collective modes that are derived using a different method.

An interesting situation arises in the absence of a magnetic field (i.e., $`\gamma =0`$), where condition (26) can be met if

$`\varphi _1+\varphi _12\varphi _0=\pi .`$ (34)
This result is in accordance with a rule pointed out in Ref. , that is, when the interaction is attractive, the relative phase coherence, as implied by the constraint (34), will spontaneously emerge if more than two BECs coexist. In the present case, the attractive force applies among three spin components, as implied by the condition $`c_1>0`$ (see Eq. (9)).

When the relation (34) holds, Eq. (26) implies that $`|\zeta _0|`$ is arbitrary, so we have $`|\zeta _1|=|\zeta _1|,`$ $`|\zeta _0|^2=12|\zeta _1|^2`$. This, combined with Eq. (34), leads to a vectorial order parameter as $`\zeta =(\zeta _1,\zeta _0,\zeta _1)=e^{i\varphi _0}(e^{i(\varphi _0\varphi _1)}\mathrm{sin}\beta /\sqrt{2},\mathrm{cos}\beta ,e^{i(\varphi _0\varphi _1)}\mathrm{sin}\beta /\sqrt{2})`$. This result agrees with Eq. (5) of Ref.  and implies that the spin components can change under spatial rotation. It is surprising that the many-body ground state predicts otherwise: from Eq. (12), we find that at zero magnetic field (hence $`F=0`$) spin components cannot change under spatial rotation but must always be the same  and is given by $`\widehat{n}_{\mathrm{𝟎}\alpha }=\left(N_{𝐤\alpha }^{}\widehat{n}_{𝐤\alpha }\right)/3`$ for $`\alpha =0,\pm 1`$, where we take into account the depletion of the condensate $`_{𝐤\alpha }^{}\widehat{n}_{𝐤\alpha }`$ due to interactions. This is a consequence of the fact that the true ground state is composed of spin-singlet “pairs” and is therefore invariant under rotation. The above result for $`\widehat{n}_{\mathrm{𝟎}\alpha }`$ and the relative phase relation (34) may be used to drastically simplify the Hamiltonian (8), giving

$`\widehat{H}_{\gamma =0}^{\mathrm{AF}}={\displaystyle \frac{c_0}{2V}}N(N1){\displaystyle \frac{c_1N}{2V}}`$ (35)
$`+{\displaystyle \underset{𝐤0}{}}[(ϵ_𝐤+c_0n)\widehat{d}_𝐤^{}\widehat{d}_𝐤+(ϵ_𝐤+c_1n)(\widehat{s}_𝐤^{}\widehat{s}_𝐤+\widehat{q}_𝐤^{}\widehat{q}_𝐤)`$ (36)
$`+{\displaystyle \frac{n}{2}}(c_0\widehat{d}_𝐤\widehat{d}_𝐤+c_1\widehat{s}_𝐤\widehat{s}_𝐤+c_1\widehat{q}_𝐤\widehat{q}_𝐤+\mathrm{h}.\mathrm{c}.)],`$ (37)
where

$`\widehat{d}_𝐤`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{3}}}(\widehat{a}_{𝐤1}e^{i\varphi _1}+\widehat{a}_{𝐤0}e^{i\varphi _0}+\widehat{a}_{𝐤1}e^{i\varphi _1}),`$ (38)
$`\widehat{s}_𝐤`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}(\widehat{a}_{𝐤1}e^{i\varphi _1}\widehat{a}_{𝐤1}e^{i\varphi _1}),`$ (39)
$`\widehat{q}_𝐤`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{6}}}(\widehat{a}_{𝐤1}e^{i\varphi _1}2\widehat{a}_{𝐤0}e^{i\varphi _0}+\widehat{a}_{𝐤1}e^{i\varphi _1}).`$ (40)
We thus find that the density, spin, and quadrupolar fluctuations provide independent excitations. The novelty of these modes is that they are phase-locked to the spin components of the condensate. Diagonalizing the Hamiltonian (37), we obtain the following excitation spectra:

$`E_{d,𝐤}^{\mathrm{AF}}`$ $`=`$ $`\sqrt{ϵ_\mathrm{k}[ϵ_\mathrm{k}+2c_0n]},`$ (41)
$`E_{s,𝐤}^{\mathrm{AF}}`$ $`=`$ $`E_{q,𝐤}^{\mathrm{AF}}=\sqrt{ϵ_\mathrm{k}[ϵ_\mathrm{k}+2c_1n]}.`$ (42)
In the thermodynamic limit, these collective modes survive only at exactly zero magnetic field, and may be viewed as singular. However, in mesoscopic situations, e.g., when the system is confined in a quasi-one-dimensional torus (to which the present theory applies with minor modifications), these modes survive at small magnetic fields, provided the $`𝐤=\mathrm{𝟎}`$, $`\alpha =0`$ mode is macroscopically occupied (for a more precise definition of “small” magnetic fields, see Ref. ).

Recently, an MIT group has observed the formation of metastable spin domains . They first prepared all atoms in the $`\alpha =1`$ state and then placed half of them in the $`\alpha =0`$ state by irradiating the rf field. Letting the system evolve freely while using the quadratic Zeeman effect to prevent the $`\alpha =1`$ component from appearing, they found that spin domains formed with the two components alternatively aligned. It was discussed  that this phenomenon is due to the imaginary frequencies of the excitation modes. Here we use our theory to confirm this hypothesis and to derive a general expression for the dispersion relation. The experimental conditions in effect amount to setting

$`|\zeta _1|^2=|\zeta _0|^2={\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{2N}}{\displaystyle \underset{𝐤\mathrm{𝟎}}{}}(\widehat{n}_{𝐤1}+\widehat{n}_{𝐤0}),\zeta _1=0.`$ (43)
Our Hamiltonian (8) reduces to

$`\widehat{H}={\displaystyle \frac{c_0+c_1}{2}}n(N1){\displaystyle \frac{c_1}{8}}nN`$ (44)
$`+{\displaystyle \underset{𝐤\mathrm{𝟎}}{}}[(ϵ_𝐤+{\displaystyle \frac{c_0n}{2}}+{\displaystyle \frac{3c_1n}{4}})\widehat{n}_{𝐤1}+(ϵ_𝐤+{\displaystyle \frac{c_0n}{2}}{\displaystyle \frac{c_1n}{4}})\widehat{n}_{𝐤1}`$ (45)
$`+{\displaystyle \frac{c_0}{8}}n(\widehat{a}_{𝐤0}\widehat{a}_{𝐤0}+\mathrm{h}.\mathrm{c}.)`$ (46)
$`+{\displaystyle \frac{c_0+c_1}{8}}n(\widehat{a}_{𝐤1}\widehat{a}_{𝐤1}+\widehat{a}_{𝐤1}\widehat{a}_{𝐤0}+\widehat{a}_{𝐤1}^{}\widehat{a}_{𝐤0}+\mathrm{h}.\mathrm{c}.)].`$ (47)
Diagonalizing the Hamiltonian, we find the dispersion relations as

$`(E_𝐤)^2`$ $`=`$ $`ϵ_𝐤^2+(2u+v)ϵ_𝐤+{\displaystyle \frac{3}{4}}v^2\pm [4(u^2+2uv+2v^2)ϵ_𝐤^2`$ (49)
$`+2v^2(2u+v)ϵ_𝐤v^3(u+{\displaystyle \frac{3}{4}}v)]^{\frac{1}{2}},`$
where $`uc_0n/2`$ and $`vc_1n/2`$. For parameters of the MIT experiment $`n10^{14}`$ cm-3, $`a_\mathrm{t}29`$Å, and $`a_\mathrm{s}26`$Å , we find that $`v/u1`$. Ignoring in Eq. (49) higher-order powers of $`v/u`$, we obtain

$`(E_𝐤)^2=ϵ_𝐤^2+(2u+v)ϵ_𝐤\pm 2(u+v)ϵ_𝐤,`$
where the plus sign corresponds to the density wave, while the minus sign corresponds to the spin wave. We see that the energy of the spin wave becomes pure imaginary for $`ϵ_𝐤<v`$, implying the formation of spin domains. The corresponding wavelength defines the characteristic length scale of the spin domains as

$`\lambda _\mathrm{c}=\sqrt{{\displaystyle \frac{3\pi }{(a_ta_s)n}}}.`$ (50)
This result agrees with that of Ref.  except for a numerical factor. Using the above parameters, we obtain $`\lambda 18\mu m`$, in reasonable agreement with the observed value of about 40$`\mu m`$ .

In conclusion, we have presented a versatile many-body theory of spin-1 dilute bose gas. Depending on the parameters, a single effective Hamiltonian (8) describes various many-body ground states, metastable spin domains, and the corresponding excitation spectra. An extension to spin-2 BEC can be similarly carried out and will be reported elsewhere.

This work was supported by a grant for Core Research for Evolutional Science and Technology (CREST) of the Japan Science and Technology Corporation (JST), by a Grant-in-Aid for Scientific Research (Grant No. 11216204) by the Ministry of Education, Science, Sports, and Culture of Japan, and by the Toray Science Foundation.


text:
# Calculation of thermodynamic properties of finite Bose-Einstein systems

## Abstract

We derive an exact recursion formula for the calculation of thermodynamic functions of finite systems obeying Bose-Einstein statistics. The formula is applicable for canonical systems where the particles can be treated as noninteracting in some approximation, e.g. like Bose-Einstein condensates in magnetic traps. The numerical effort of our computation scheme grows only linear with the number of particles.
As an example we calculate the relative ground state fluctuations and specific heats for ideal Bose gases with a finite numbers of particles enclosed in containers of different shapes.

With the observation of Bose-Einstein condensation (BEC) of magnetically  and optically  trapped atoms new insights into the nature of this state of matter have been given. The experimental situation is in all cases quite different from the ideal gas treated within the grand canonical ensemble, which is the standard textbook example. First, the number of particles within the traps is fixed and finite, which suggests a canonical or microcanonical treatment of the systems. Second, the confining trap potentials greatly influence the condensate properties. Third, although the trapped gases are quite dilute the validity of the treatment as non-interacting particle gases has to be checked from case to case.
Even within the approximation of non-interacting particles the calculation of the thermodynamic properties of the Bose-Einstein systems remains a hard to tackle mathematical problem. Recently, some approximate methods to calculate the fluctuation of the ground state occupation number in a trapped Bose-Einstein condensate have been developed . Here we present an exact method to calculate all thermodynamic quantities of finite canonical Bose systems, given the one particle density of states.

As the starting point we utilize the recursive formula of the canonical partition function for a system of $`N`$ noninteracting bosons as given in :

$$Z_N(\beta )=\frac{1}{N}\underset{k=1}{\overset{N}{}}Q_k(\beta )Z_{Nk}(\beta ),$$
(1)
where $`Q_k(\beta )=Z_1(k\beta )=_i\mathrm{exp}(k\beta ϵ_i)`$ is the one-particle partition function at the temperatures $`k\beta `$ and $`Z_0(\beta )=1`$. The microcanonical partition $`\mathrm{\Gamma }_N(E)`$ can be calculated by an inverse Laplace-transform of (1) and is given by

$`\mathrm{\Gamma }_N(E)`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{k=1}{\overset{N}{}}}{\displaystyle \frac{1}{2\pi i}}{\displaystyle _{ci\mathrm{}}^{c+i\mathrm{}}}d\beta \mathrm{exp}(\beta E)Q_k(\beta )Z_{Nk}(\beta )`$ (2)
$`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{k=1}{\overset{N}{}}}{\displaystyle _0^E}dE^{}\mathrm{\Gamma }_1^k(E^{})\mathrm{\Gamma }_{Nk}(EE^{}),`$ (3)
where $`\mathrm{\Gamma }_1^k(E)`$ is the inverse Laplace-transform of $`Q_k(\beta )`$ and $`\mathrm{\Gamma }_0(E)=\delta (E)`$. A similar, slightly less general, equation has recently been derived by Weiss and Wilkens .

Eq.(1) can be used to calculate all thermodynamic quantities by appropriate differentiation of $`\mathrm{ln}Z_N`$. However, in any case $`Z_N`$ occurs as a normalization factor and has to be calculated explicitly. This turns out to be a major drawback. First, the numerical effort to calculate $`Z_N`$ grows with the square of the particle number $`N`$. Moreover, since $`Z_N(\beta )`$ grows exponentially with $`N`$ multiple precision arithmetic is required for proper calculation. We will present a method avoiding these difficulties.

To ease our derivations we rewrite $`Z_N(\beta )`$ utilizing the $`𝒵`$-transform and define:

$`𝒵(Z)`$ $`=`$ $`F(x)={\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{Z_k(\beta )}{x^k}},`$ (4)
$`𝒵(Q)`$ $`=`$ $`G(x)={\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{Q_k(\beta )}{x^k}},`$ (5)
where we define $`Q_0(\beta )=0`$. Taking advantage of the basic properties of the $`𝒵`$-transform eq.(1) can be written in the form

$$x\frac{\mathrm{d}}{\mathrm{dx}}F(x)=F(x)G(x),$$
(6)
yielding ***Note that $`F(x)`$ is closely related to the grand canonical partition function.

$$F(x)=\mathrm{exp}\left(\underset{k=1}{\overset{\mathrm{}}{}}\frac{Q_k(\beta )}{k}x^k\right).$$
(7)
Applying the inverse $`𝒵`$-transform we may write $`Z_N(\beta )`$ as

$$Z_N(\beta )=\frac{1}{2\pi i}_CF(x)x^{N1}dx,$$
(8)
where $`C:=\{zC:x=r\}`$ and $`r`$ has to satisfy the condition $`Z_N(\beta )\mathrm{exp}(rN)`$. Alternatively we may write:

$$Z_N(\beta )=\frac{1}{N!}\frac{\mathrm{d}^\mathrm{N}}{\mathrm{dx}^\mathrm{N}}F(1/x)|_{x=0}.$$
(9)
Using eq.(9) the number of particles with energy $`ϵ_i`$ can be calculated by

$`\eta _i(N,\beta )`$ $`=`$ $`{\displaystyle \frac{1}{\beta }}{\displaystyle \frac{}{ϵ_i}}\mathrm{ln}Z_N(\beta )`$ (10)
$`=`$ $`{\displaystyle \frac{1}{Z_N(\beta )}}{\displaystyle \underset{k=1}{\overset{N}{}}}\mathrm{exp}(\beta kϵ_i)Z_{Nk}(\beta )`$ (11)
Some reordering yields:

$$\eta _i(N+1,\beta )=\frac{Z_N(\beta )}{Z_{N+1}(\beta )}\mathrm{exp}(\beta ϵ_i)(\eta _i(N,\beta )+1).$$
(12)
Since the particle number is a conserved quantity in the canonical ensemble the direct calculation of the normalization factor can be omitted by using the relation

$$\frac{Z_N(\beta )}{Z_{N+1}(\beta )}=\frac{N+1}{_{i=0}^{\mathrm{}}\mathrm{exp}(\beta ϵ_i)(\eta _i(N,\beta )+1)}.$$
(13)
For Fermi systems the following recursion formula

$$\eta _i(N+1,\beta )=\frac{Z_N(\beta )}{Z_{N+1}(\beta )}\mathrm{exp}(\beta ϵ_i)(1\eta _i(N,\beta )),$$
(14)
with

$$\frac{Z_N(\beta )}{Z_{N+1}(\beta )}=\frac{N+1}{_{i=0}^{\mathrm{}}\mathrm{exp}(\beta ϵ_i)(1\eta _i(N,\beta ))}.$$
(15)
can be derived in a similar manner. In practice, only a limited number of energy levels has to be taken into account, since the occupation probability rapidly decreases with increasing energy eigenvalues. Equations (12) and (13) are extremely useful in practical calculations. The numerical effort to calculate the occupation numbers grows only linear with the number of particles. Moreover, only a moderate arithmetic precision is required. Having the occupation probabilities at hand the energy expectation value is given by:

$$E(N,\beta )=\underset{i=0}{\overset{\mathrm{}}{}}ϵ_i\eta _i(N,\beta ).$$
(16)
The calculation of the fluctuation of the occupation probabilities $`\delta \eta _i(N,\beta )`$ is a little bit more complicated and contains another recursion:

$`\left(\delta \eta _i(N+1,\beta )\right)^2={\displaystyle \frac{1}{\beta ^2}}{\displaystyle \frac{^2}{ϵ_i^2}}\mathrm{ln}(Z_{N+1}(\beta ))`$ (17)
$`=`$ $`{\displaystyle \frac{1}{\beta }}{\displaystyle \frac{}{ϵ_i}}\eta _i(N+1,\beta )`$ (18)
$`=`$ $`{\displaystyle \frac{Z_N(\beta )}{Z_{N+1}(\beta )}}\mathrm{exp}(\beta ϵ_i)\{\delta ^2\eta _i(N,\beta )`$ (19)
$`+`$ $`[\eta _i(N+1,\beta )+1][\eta _i(N,\beta )\eta _i(N+1,\beta )+1]\}.`$ (20)
To illustrate the usefulness of our recursion formulas we consider the ideal gas with parameters of liquid Helium in containers of different shapes :
i) a cube with side-length $`L_x,L_y,L_z`$, and energy levels

$$E_{n_x,n_y,n_z}=\frac{\pi ^2\mathrm{}^2}{2m_{He}}\left(\frac{n_x^2}{L_x^2}+\frac{n_y^2}{L_y^2}+\frac{n_z^2}{L_z^2}\right)$$
(21)
ii) a sphere with radius $`a`$ and energy levels

$$E_{n,l}=\frac{\mathrm{}^2}{2m_{He}a^2}u_{n,l}$$
(22)
and degeneracy $`\sigma _{n,l}=2l+1`$, and
iii) a cylinder with diameter $`d=2a`$, height $`L`$, energy levels

$$E_{n,l,m}=\frac{\mathrm{}^2}{2m_{He}}\left(\frac{v_{n,l}^2}{a^2}+\frac{m^2\pi ^2}{L^2}\right)$$
(23)
with $`n=1,2,\mathrm{}`$, $`l=1,2,3,\mathrm{}`$, and $`m=\mathrm{},1,0,1,\mathrm{}`$.
We denoted the zeros of the half integer Bessel functions $`J_{n+1/2}(r)`$ by $`u_{n,l}`$ and the zeros of the integer Bessel function $`J_n(r)`$ by $`v_{n,l}`$.

Fig. 1 displays the specific heats and the fluctuations of the ground state occupation number $`\delta \eta _0/N`$ as a function of the canonical temperature for different trap geometries and $`N=100,1000`$, and $`10000`$ He-atoms. In all cases the particle density is taken to be $`\rho =0.0216\AA ^3`$. With growing system size the differences between the specific heats for the different trap geometries almost vanish and approach the typical shape of the curve for the ideal Bose gas. I.e., with respect to the specific heat the boundary conditions get more and more unimportant with increasing volume. In contrast, the ground state fluctuations exhibit a completely different behaviour. The cubic box, the compact cylinder with equal diameter and height, and the sphere show almost equal ground state fluctuations for all system sizes, while the ground state fluctuations of the stretched box and the stretched cylinder are remarkably larger for temperatures below the critical temperature. This effect is not unexpected because restricting the particle motion in one or two dimensions makes the system act like a lower dimensional system, which are known to show larger fluctuations. Moreover, the differences between the fluctuations of the stretched traps and the compact traps do not decrease with increasing system size. The reason for this behaviour is found in the energy difference between the ground state and the first excited level, which is much larger for the stretched traps than for the compact traps. Since $`\delta \eta _0/N`$ decreases approximately with $`N^{1/3}`$ the infinite particle number limit is the same for all trap geometries. Under experimental considerations our results imply that the stability of the condensate fraction in anisotropic traps should be considerably smaller than in isotropic traps. In Fig. 1 we plotted $`\delta \eta _0/N`$ to allow good comparison with previous published results . Since this quantity goes to zero as the system size increases it is a bad indicator for phase transitions. The relative ground state fluctuation $`\delta \eta _0/\eta _0`$ shown in Fig. 2 is much more conclusive in this respect.

True phase transition only occur for infinite systems. However, it is well known from other systems, e.g. finite spin lattices and clusters , that finite systems already display the onset of phase transitions. Instead of having a well defined critical temperature the transition occurs in a broader crossover region. As can been extracted from Fig. 2 the crossover region, which is indicated by the sharp increase of $`\delta \eta _0/\eta _0`$, extends even for the cubic box with 10000 atoms over a temperature range of 0.5 K. For the stretched trap the crossover region is about twice as large.

All calculations have been performed on an IBM-43P (233 MHz) workstation. For 1000 particles a run with 200 temperature points took about six minutes, for 10000 particles one hour. A calculation of the same quantities with the recursion formula given in equation (1) takes for 1000 particles about six hours (and would take at least 600 hours for 10000 particles). Moreover, due to the numerical instabilities connected with (1) $`N=1000`$ was the largest particle number we achieved with this formula, even though we utilized a multiple precision package.

We expect the recursion formulas to be quite useful for calculating the properties of dilute atomic Bose gases in magnetic traps with different geometry. For this propose we provide an easy to use Java programThe code is available within the World-Wide-Web:
http:$`\backslash \backslash `$www.physik.uni-oldenburg.de$`\backslash `$borrmann$`\backslash `$BEC, which requires as input only the energy level distribution with appropriate degeneracy and calculates the basic properties of finite Bose systems with particle numbers up to 10000. A slightly faster Fortran code is available upon Email requestCorresponding author: borrmann@uni-oldenburg.de.