text:
# Nucleation and Growth of GaN/AlN Quantum Dots

## I Introduction

Zero-dimensional semiconductor quantum dots (QDs) have attracted much interest in the last decade due to their multiple potential applications ranging from low-threshold lasersArakawa ; Bimberg  via single-electron tunneling devicesOdintsov ; Kouwenhoven  to possible realizations of qubits for quantum computation.QC1 ; QC2  A versatile method for the fabrication of semiconductor QDs is their self-assembled growth following the Stranski-Krastanov (SK) growth mode.ReviewsSK  This mode usually occurs during the growth of semiconductor epilayers under compressive strain. Example material systems are InxGa1-xAs/GaAs,Goldstein ; Leonard93 ; Moison  SixGe1-x/Si,Eaglesham ; Mo  CdSe/ZnSe,Xin  or GaN/AlN.Daudin ; Widmann  In this mode, atoms are initially deposited in form of a two-dimensional pseudomorphic wetting layer. The associated strain energy increases with the thickness of the wetting layer and is finally elastically relieved by the formation of islands.

The usefulness of such self-assembled nanostructures relies on the ability to obtain homogeneous size distributions as well as to control their size, density and position. Many theoretical contributions have enhanced our understanding of the size distributions of SK-grown islands, but some controversy remains.Shchukin ; Daruka ; Williams ; Priester ; Chen ; Ross ; Jesson  Also, despite many experimental studies, the influence of growth parameters on the size and density of such islands is not fully understood, owing to the complexity of the physics of strained layer growth.

In this work, we present results on the nucleation of GaN islands on AlN following an SK mode, in particular on the dependence on the amount of deposited GaN and the substrate temperature. We find that the qualitative behavior is similar to that found in other systems, e.g. for InAs/GaAs and Ge/Si. This reinforces the idea that there are common features of semiconductor SK growth, which are rather universal and independent of the specific material system. Absolute island sizes then densities will depend on the material system, possibly through material parameters like lattice misfit, elastic constants, or surface energies. Hence, we may further discuss the experimental data in the framework of an equilibrium model based on statistical physics that was originally developed in work on the Ge/Si system.Rudd

## II Experimental

The samples have been grown in a MECA2000 molecular-beam epitaxy (MBE) chamber equipped with standard effusion cells for Ga and Al evaporation. The chamber also contains an rf plasma cell to provide active nitrogen for GaN and AlN growth. The pseudo-substrates used were about 2 $`\mu `$m thick GaN(0001) (Ga-polarity) layers grown by metal-organic chemical vapor deposition on sapphire. The substrate temperature $`T_S`$ was measured by a thermocouple in mechanical contact to the backside of the molybdenum sample holder. To ensure substrate temperature reproducibility, each series of samples described below was grown on a single molybdenum substrate holder.

Prior to all experiments, a 100 nm thick GaN layer was grown under Ga-rich conditions on the pseudosubstrates to avoid the influence of a possible surface contamination layer. Subsequently, a 300 nm thick AlN film was deposited under Al-rich conditions at a substrate temperature of 730 C. We have found by reflection high-energy electron diffraction (RHEED) and high-resolution X-Ray diffraction that this thickness is sufficient for the AlN layer to be virtually fully relaxed with a residual in-plane strain of $`ϵ_1<0.1`$%.Bellet  An atomic-force microscopy (AFM) image of such an AlN(0001) surface is shown in Fig. 1. The surface is characterized by about 30 nm wide terraces and spiral hillocks, similar to GaN surfaces grown under equivalent conditions.AdelmannJAP

The growth rate and the GaN coverage have been experimentally determined for each sample by RHEED oscillations occurring during the growth of the wetting layer prior to island formation.Mula  Typically, growth rates for different layers (and different samples) were reproducible within about 0.01 monolayers (ML)/s. The GaN coverage was then calculated with a precision better than 0.1 ML. It is worth noting that no measurable influence of the growth rate on GaN island properties has been found in the range between 0.1 and 0.6 ML/s. As, on the other hand, the Ga/N flux ratio has been found to have a crucial influence on the growth mode,Mula  it was fixed to 0.8 (N-rich conditions), which leads to a critical thickness for the SK transition of 2.25 ML.Mula ; AdelmannThesis  The critical thickness was measured for each sample and found to be reproducible within 0.1 ML.

To study the GaN islands *as grown*, the samples were rapidly quenched to room temperature under an N-flux. *Ex-situ* AFM was then used to study the GaN morphology after exposure of the samples to air.

## III Results

### III.1 Dependence on GaN coverage

The parameter most directly governing the properties of GaN islands grown in an SK mode on AlN is the GaN coverage $`\mathrm{\Theta }`$, which effectively describes the time evolution of the islands during growth. To study this evolution, a series of samples has been grown at a substrate temperature of 730 C and a growth rate of 0.15 ML/s. The GaN coverage was varied between $`\mathrm{\Theta }=1.8`$ ML and $`\mathrm{\Theta }=4.6`$ ML.

Figure 2(a) shows an AFM image of the morphology of GaN layers obtained after the deposition of 1.8 ML, *i.e.* for a coverage well below the critical thickness of 2.25 ML. We find that the morphology is unchanged with respect to that of the AlN pseudosubstrate: the surface is characterized by about 30 nm wide terraces and spiral hillocks. We can thus infer that the growth of about the first 2 ML of GaN occurs in a layer-by-layer mode since RHEED oscillations are observed.Mula

When the second monolayer is completed, the morphology changes, as evidenced in Fig. 2(b) for $`\mathrm{\Theta }=2.2`$ ML, *i.e.* immediately before the SK transition. Remainders of terraces and spiral hillocks are still visible, but the surface is characterized on a short scale by 1–2 ML high flat 2D islands with typical diameters of $`d20`$ nm (see Fig. 2(b)). The behavior can thus be described by a transition from a layer-by-layer growth to multilayer growth at around 2.0 ML. Such an occurrence of 2D precursor islands prior to the genuine 2D-3D SK transition has also been observed in the InAs/(Al,Ga)As system.Cirlin ; Kitabayashi ; Ramachandran ; Ballet

When the growth is continued, the genuine 2D-3D transition occurs and 3D GaN islands are formed as shown in Fig. 3. We find that their density increases strongly with GaN coverage and saturates around 3.0 ML at a value of $`1.8\times 10^{11}`$ cm-2 \[see Figs. 3(a) and (b)\]. Further GaN deposition does not lead to an increase in island density but instead islands grow in size. However, the islands do not grow continuously in size but a bimodal size distribution is observed \[see Figs. 3(c) and (d)\].

The variation of the island size is summarized in Figs. 4 and 5. We see that, initially, the islands’ height increases from 0.7 nm to 1.6 nm between 2.4 and 2.8 ML of GaN coverage. For higher GaN coverage, a bimodal size distribution is observed with the average height of the first mode islands remaining constant at 1.6 nm independent of coverage. In contrast, the average height corresponding to the second mode increases continuously reaching 4.2 nm at 4.7 ML coverage with no sign of saturation in the examined coverage range. Differences in the shapes of mode 1 and mode 2 islands will be discussed below.

The variation of the island density is depicted in Fig. 6. We observe that the total density increases strongly after the 2D–3D transition but saturates after 2.8 ML, i.e. at the coverage where the island size distribution becomes bimodal. The partial densities of the two modes are shown in the inset in Fig. 6. For coverages below 2.8 ML, the density of mode 1 islands is identical with the total density and increases strongly with coverage. After the deposition of 2.8 ML, the density of mode 2 islands increases strongly, similarly to the behavior of mode 1 islands after the 2D–3D transition, whilst the density of mode 1 islands decreases. As the total island density remains approximately constant, mode 1 islands transform into mode 2 islands, probably without additional nucleation of new (mode 1) islands.

### III.2 Dependence on substrate temperature

The influence of the substrate temperature is studied in a series of samples with $`\mathrm{\Theta }=3.0`$ ML deposited at substrate temperatures between $`T_S=690`$C and $`T_S=760`$C. At lower temperatures, no SK growth is observed and GaN grows in a pseudo-2D mode,Daudin ; Mula ; AdelmannThesis  probably due to excessively low adatom mobility.Zywietz  Higher substrate temperatures are prohibited by our experimental apparatus due to indium bonding of the substrates to the substrate holder. All samples in this series were grown at a growth rate of 0.25 ML/s.

An AFM image of the morphology of two GaN layers grown at substrate temperatures of $`T_S=700`$C and $`T_S=760`$C, respectively, are shown in Fig. 7. We observe that the total island density decreases rapidly with substrate temperature, in keeping with previous resultsWidmann  and also with results obtained for the InAs/GaAsLeon1  and Ge/SiKamins  systems. The variation of the GaN island density as a function of substrate temperature for a coverage of $`\mathrm{\Theta }=3.0`$ ML is depicted in Fig. 8. We find an approximately exponential decrease of the island density with increasing substrate temperature from $`4.1\times 10^{11}`$ cm-2 at $`T_S=690`$C to $`8.0\times 10^{10}`$ cm-2 at $`T_S=760`$C.

As we have seen above, the island density tends to saturate at sufficiently high GaN coverage. This saturation density is shown as a function of substrate temperature in the inset in Fig. 8. For substrate temperatures of $`T_S730`$C, the saturation density is similar to the density for $`\mathrm{\Theta }=3.0`$ ML, in agreement with the results in the last section. However, for higher substrate temperature, the saturation density is significantly higher, indicating that saturation occurs for GaN coverages larger than 3.0 ML.

Figure 9 shows the island height distribution for $`\mathrm{\Theta }=3.0`$ ML and substrate temperatures as indicated. The behavior is summarized in Fig. 10. We find a single approximately Gaussian distribution (in the limit of the statistical precision) at a height of 1.6 nm for substrate temperatures $`T_S720`$C. These islands correspond to the mode 1 islands observed at $`T_S=730`$C and discussed above. At low temperature, the distribution is thus still monomodal after 3.0 ML. For $`T_S=730`$C, a shoulder appears at the high island side of the distribution and transforms into a clearly separated second mode at higher temperatures. Hence, we find that bimodal distributions occur earlier (for lower GaN coverage) at higher substrate temperature. Another remarkable finding is that the height of the mode 1 islands remains constant as a function of temperature, whereas the height of mode 2 islands increases with temperature, again without signs of saturation in the examined temperature range.

### III.3 Discussion

The above results demonstrate the occurrence of bimodal island size distributions at high GaN coverage and/or substrate temperature; at low GaN coverage and/or substrate temperature, the size distribution is monomodal. Remarkably, when bimodal distributions occur, the size of mode 1 islands (smaller islands) appears independent of growth parameters such as GaN coverage and substrate temperature.

As discussed above, the variation of partial densities of both modes implies that first mode 1 islands are nucleated, their density saturates and the growth proceeds further by transformation of mode 1 islands into mode 2 islands. Thus the question arises, whether this transformation involves a shape change of the islands similar to the behavior observed in the Ge/Si system, where pyramidal islands are observed to transform into “dome”-shaped islands during growth.Medeiros ; Williams  Some information about the island shape can be gathered by RHEED. The RHEED images corresponding to pure mode 1 island morphologies and morphologies with a mixture of mode 1 and mode 2 islands (but with a majority of mode 2 islands) are both characterized by the same $`\left\{10\overline{1}3\right\}`$ facets, showing a six-fold rotation symmetry. It thus appears that both types of islands are (truncated) hexagonal pyramids with $`\left\{10\overline{1}3\right\}`$ sidewalls, in agreement with previous results.Daudin ; Widmann  However, it is clear that the shape of such islands cannot be obtained from the RHEED pattern alone and further work is necessary to determine the precise shape of both types of islands.

A more detailed analysis can be done by extracting the aspect ratio $`\tau `$ of the islands from the AFM data. Figure 11(a) depicts the aspect ratio of mode 1 islands for two different GaN coverages (at which monomodal height distributions are observed) at a substrate temperature of $`730`$C. In spite of the absolute increase of the islands’ size, their aspect ratio remains constant and we find $`\tau =0.12`$. As shown in Fig. 11(b), the aspect ratio of mode 1 islands remains constant ($`\tau =0.13`$) after the transition to a bimodal distribution but the mode 2 islands have a significantly larger aspect ratio of $`\tau =0.19`$. Both aspect ratios correspond to truncated pyramids remark  but mode 1 islands are flatter.

From the aspect ratio data, we can further calculate the amount of GaN contained in the 3D islands as a function of deposited GaN. The result is shown in Fig. 12. We find that the data are consistent with a 2 ML thick wetting layer, independent of the amount of deposited GaN. It thus appears that in this system, the wetting layer is stable and does not contribute to island growth. The deviation for GaN coverages just after the critical thickness of 2.25 ML can be explained by the presence of 2D precursor islands, which still contain significant material but are not taken into account as 3D islands.

## IV Modeling

In order to understand more precisely why the observed island distributions form, we analyze the nitride nanostructure growth in terms of a thermodynamic equilibrium model that has so far only been applied extensively to the Ge/Si system.Rudd  To date, little effort has been made to fit data on size distributions of nitride nanostructures to any model. Nakajima *et al.*Nakajima  have undertaken a purely theoretical study and explicitly derived thickness-composition phase diagrams for the expected growth mode for GaInN/GaN and GaInN/AlN, using known materials parameters. However they did not take into account the shapes and sizes of experimentally observed nitride nanostructures, nor did they seek to predict size distributions. Here we take a rather different approach. We apply the thermodynamic model to the AFM data on nitride nanostructure distributions grown at one particular temperature in order to assess whether the bimodal distributions of nanostructures are consistent with thermodynamic equilibrium, as opposed to being configurations along an unstable ripening trajectory. The experimental distributions are described reasonably well by the model, as explained below. We then apply the model to understand the variations in the bimodal distributions with growth temperature at fixed coverage, and again find reasonable agreement over a range of growth temperatures.

The thermodynamic equilibrium model makes explicit predictions for the form of the bimodal epitaxial nanostructure distributions. The model is governed by a set of parameters that describes the size dependence of the internal energy of an epitaxial nanostructure. Following Shchukin *et al.*,Shchukin  we describe the internal energy $`ϵ_i`$ of the $`i`$th individual island of type $`X`$ and volume $`\nu _i`$ as follows:

$$ϵ_i=A_X\nu _i+B_X\nu _i^{2/3}+C_X\nu _i^{1/3}+D_X.$$
(1)
Using this expression for the internal energy, the island size distribution $`f(\nu _i)`$ is given by

$$f(\nu _i)=\mathrm{exp}\left(\frac{ϵ_i\mu \nu _i}{k_BT_S}\right),$$
(2)
where $`\mu `$ is the chemical potential and $`k_B`$ is Boltzmann’s constant.

Equation (1) describes the chemical and elastic contributions to the island energy. The first term arises from bulk strain, the second from surface and interfacial energies and island-island interactions, and the third from surface stress and edge energies (suppressing the log dependence of the edge elastic relaxation energy). The equation also includes a size-independent term that is relevant only insofar as it differs from one nanostructure type to another. For the internal energy (1), a minimum may exist if $`B`$ and $`C`$ have opposite signs. There has been some debate over which of the two constants is likely to be positive and which negative. In the initial development of a “shape-map” for the bimodal Ge/Si system by Williams *et al.*,Williams  $`B`$ was assumed to be negative and $`C`$ positive. In more recent work on Ge/Si, Rudd *et al.*Rudd  have suggested that $`B`$ should be positive, since a Ge(001) surface does not spontaneously roughen, though they also acknowledge that, since 3D growth is preceded by the appearance of increasing numbers of defects in the wetting layer with their own associated energy, the behaviour of the Ge(001) surface of a bulk crystal may not be relevant and $`B`$ could be negative. Here, we find that the option most consistent with the nanostructure size distributions observed by AFM is $`B>0`$ and $`C<0`$.

Unlike in the most recent developments of the equilibrium model,Rudd  we have not, as yet, explicitly included the effects of the elastic interaction between pairs of islands. We also have not undertaken a self-consistent calculation of the chemical potential. Instead, we have taken a similar approach to Williams *et al.*:Williams  Using the expressions for the distributions (2) and internal energies (1), we have fitted the model to the available data at several coverages and temperatures, in order to determine approximately how the constants $`A^{}`$, $`B`$, $`C`$, and $`D`$ vary as a function of growth parameters. For convenience, we have introduced $`A^{}=A\mu `$. We then used our empirical functions to fit the model across multiple data sets, to determine how well it compares with the overall behaviour of the system. By comparison with the more sophisticated work of Rudd *et al.*Rudd  on Ge/Si, we expect $`B`$ to be dependent on the coverage and $`A^{}`$ to depend on the chemical potential, a non-trivial function of coverage and temperature.

### IV.1 Evolution of size distribution with coverage

In investigating the variation of the parameters $`A^{}`$, $`B`$, $`C`$, and $`D`$ as a function of coverage, we found that $`C`$ and $`D`$ could be treated as being independent of coverage, but that $`B`$ increased monotonically with increasing coverage. This coverage dependence is consistent with elastic island-island interactions, which are expected to scale with $`\nu _i^{2/3}`$, and which will increase with increasing coverage. $`B\left(\mathrm{\Theta }\right)`$ could be approximated to a straight line with non-zero slope for both the smaller ($`S`$) and the larger ($`L`$) island types. Thus (for $`X\{S,L\}`$):

$$B_X(\mathrm{\Theta })B_{0X}+b_X\mathrm{\Theta }.$$
(3)
Similarly, Williams *et al.*Williams  found that, for the Ge/Si case, $`C`$ was independent of coverage and $`B`$ varied linearly with coverage. However, they were also able to approximate the variation of $`A^{}`$ with coverage as a straight line with non-zero slope. This does not appear to be valid in this case, and we have not found a simple functional relationship between $`A^{}`$ and coverage. Hence in fitting the model across multiple data sets we have calculated a separate value of $`A^{}`$ for each data set, but have assumed a constant relationship between the value of $`A^{}`$ for the larger and smaller islands. Essentially:

$$A_L^{}(\mathrm{\Theta })=g(\mathrm{\Theta })\mathrm{and}A_S^{}(\mathrm{\Theta })=g(\mathrm{\Theta })+A_0,$$
(4)
where $`A_0=A_SA_L`$, representing the intrinsic differences in bulk elastic energy between the two island types. The function $`g(\mathrm{\Theta })`$, which has not been found explicitly, represents the coverage dependence of the chemical potential.

Hence, in fitting across four coverages at which large and small islands coexist, we used the parameters $`B_{0S}`$, $`B_{0L}`$, $`b_S`$, $`b_L`$, $`C_S`$, $`C_L`$, $`D_S`$, $`D_L`$, $`A_0`$, and four separately fitted values of $`g(\mathrm{\Theta })`$. This clearly gives us a large number of parameters—however, the model must fit a large number of observable features of the data. We are able to provide a reasonable fit for each shape for the relative height, position, width and skewness of the volume distributions (essentially 15 experimental observables for each island shape). Example fits are shown in Fig. 13. The fitting has been done using a $`\chi ^2`$ approach and the goodness-of-fit may be judged by calculating an overall $`Q`$-value across the data sets, $`Q_{\mathrm{total}}`$. For the data sets under consideration, $`Q_{\mathrm{total}}=0.22`$. Since the available data sets are much smaller for this system than for the Ge/Si system (in which counting large numbers of islands is facilitated by the larger sizes of the islands, and the clear shape difference between island types) this $`Q`$-value is satisfactory. For comparison we have also attempted to fit the data with the original model of Williams *et al.*,Williams  and with some Gaussian functions. These fits were markedly worse and it was not possible to find parameters that fit the entire set of data with a reasonable $`Q`$-value. From these findings we may infer that the data are consistent with the equilibrium model; however, the model can in no way be said to have been proven. The model contains a large number of parameters for which little *a priori* information is available in the nitride systems, and the measured island distributions would benefit from more data in order to reduce the statistical uncertainty. Hence, the discussion that follows is somewhat tentative.

Using the calculated parameters we can plot the internal free energy versus volume curves for both the smaller and the larger islands. The results are shown in Fig. 14 for three different coverages. There is a significant evolution of the curves as the coverage increases. Similar results were observed by Williams *et al.*Williams  for the evolution of internal energy with coverage in the Ge/Si system, and were attributed to the interactions of the nanocrystals. The value of $`B(\mathrm{\Theta })`$ affects the width of the curve and the horizontal position of the energy minimum. In the Ge/Si case, Williams *et al.* found that for the larger islands (domes) the value of $`B(\mathrm{\Theta })`$ increases (becomes less negative) with increasing coverage, whilst the value of $`B(\mathrm{\Theta })`$ for the smaller islands (pyramids) decreases, which suggested that the repulsion between domes was stronger than the repulsion between pyramids. In contrast, the work of Rudd *et al.*Rudd  (also on Ge/Si) suggests that $`B(\mathrm{\Theta })`$ increases with increasing coverage at a similar rate for both pyramids and domes. Similarly, in the current case, we find that the $`B`$-parameter increases with increasing coverage at almost the same rate for both island types, indicating no significant difference between the elastic island-island interactions for each island type.

Changes in $`A^{}(\mathrm{\Theta })`$ tend to shift the internal-free energy versus volume curves vertically. At low coverages, the minimum internal energy is lower for the smaller islands, but at high coverages, it is lower for the larger islands. Hence the formation of large islands rather than small islands is preferred at high coverages. However, the elastic interactions between islands may mean that not all small islands grow to become large islands. Some may be destabilised by neighbouring islands and dissolve. One surprising feature of the fits is that $`A^{}`$ is found to be positive for the larger islands, for all coverages. Since the bulk elastic energy is expected to decrease when an island is formed, this suggests that the chemical potential is negative and relatively large (since $`A^{}=A\mu `$). This finding is not yet well understood.

### IV.2 Evolution of size distribution with temperature

In order to understand the variation of the size distribution with temperature, we must consider the temperature dependence of the free energy of each island type. In principle, $`A^{}`$, $`B`$, and $`C`$ may all be temperature dependent. The size-independent $`D`$-parameter is constant and is found to be the same for both small and large islands. Williams *et al.*Williams  performed some limited studies on the temperature dependence of the fitting parameters in Ge/Si and found that only $`A^{}`$ exhibited significant temperature dependence. Similarly, in this case, $`C`$ is found to be essentially independent of temperature. However, for both island types, we have found that $`B`$ exhibits a significant temperature dependence, such that for $`X\{S,L\}`$:

$$B_X(T_S)=B_{1X}+\beta _XT_S.$$
(5)
The $`B`$-parameter increases at a similar rate for the larger and the smaller islands, implying that the surface energy contribution from each island type varies in a similar way with temperature. This is congruent with the suggestion that the $`\left\{10\overline{1}3\right\}`$ facets seen in RHEED dominate both island types. These results are not, however, incompatible with the existence of other facets on the larger island types.

The $`A^{}`$-parameter, which was found to exhibit a complex dependence on coverage, also exhibits a complex dependence on temperature. Hence in fitting the model across a range of temperatures, we have calculated values of $`A^{}`$ individually for each data set, but have assumed a constant relationship between the value of $`A^{}`$ for the large islands and the value of $`A^{}`$ for the small islands, such that

$$A_L^{}(T_S)=h(T_S)\mathrm{and}A_S^{}(T_S)=h(T_S)+A_1.$$
(6)
For the larger islands, $`A^{}`$ is again found to be positive at lower temperatures but decreases smoothly with temperature and changes sign at $`755`$C. The fits are shown in Fig. 15; $`Q_{\mathrm{total}}=0.73`$. The function $`h(T_S)`$, which has not been found explicitly, includes the temperature dependence of the chemical potential, and any variation in the elastic relaxation with temperature.

This model is only valid if a quasi-equilibrium may be assumed. This assumption is particularly problematic with regards to variation in temperature, since in particular at lower temperatures, the system may not reach equilibrium, whilst at higher temperatures significant interdiffusion between the islands and the substrate may occur. This alloying will effectively reduce the mismatch between the islands and the substrate - strongly affecting the island size distributions, since the mean equilibrium nanostructure volume varies inversely as the sixth power of the mismatch.Dorsch  Rudd *et al.*’s more sophisticated work,Rudd  in creating the nanostructure diagram for the Ge/Si system, could be augmented to account for uniformly distributed alloying, and for kinetic effects. In the past, whilst the temperature dependence of the equilibrium island size distribution has been considered theoretically,Rudd ; Williams  very little high quality data has been available for comparison. The experimental studies by Williams *et al.*Williams  only involved temperatures of 550 C, 600 C, and 650 C, and there was considerable uncertainty in the temperature measurement. Additionally, the sample grown at 550 C is believed to have been significantly influenced by kinetic effects, and the sample grown at 650 C experienced significant intermixing. Hence the data considered here, relating to a smaller range of more accurately measured temperatures, are very valuable as a first quantitative test of the temperature dependence of the equilibrium model. In addition, the GaN/AlN system presents an important advantage over the Ge/Si system since it has been show that alloying is very limited even at high growth temperatures.Arlery

In considering the results of the fitting, we notice that fitting across the temperature data sets gives a higher $`Q_{\mathrm{total}}`$ than fitting across the coverage data, with an essentially similar fitting methodology. One might thus speculate that, since the samples used for the coverage data set were grown at the relatively low temperature of 730 C, full equilibrium has not been achieved and kinetic effects are having some influence. A worthwhile approach, in investigating this further, will be to examine a series of samples that have been annealed at growth temperature to allow further evolution towards equilibrium.

## V Conclusion

Using *ex situ* AFM, we have studied the evolution of thin GaN layers grown by plasma-assisted MBE on AlN(0001) layers at substrate temperatures between 700 and 750 C. Initially, 2 ML of GaN grow in a 2D mode, followed by the occurrence of 2D islands. These islands act as precursors for 3D islands, which appear after an SK transition around 2.3 ML. During further growth, in particular at higher temperatures, a bimodal island size distribution is observed. Remarkably, the size of mode 1 islands is found to be independent of coverage and temperature, whereas the size of mode 2 islands increases with coverage and temperature. The analysis of the partial island densities reveals that, whilst the total island density remains constant, mode 1 islands transform during growth into mode 2 islands. The aspect ratio of the islands is measured for both types of islands and it is found that they are characterized by distinctively different aspect ratios, whereas no additional facets are observed in the RHEED pattern due to mode 2 islands.

These findings are examined in the framework of an equilibrium model for SK growth. We find satisfactory agreement with the experimental data, suggesting that, as more data become available, extending the equilibrium approach may be helpful in understanding and tailoring nitride nanostructure distributions. We have examined the variation of each of the fitting parameters with growth conditions, and have considered how this may relate to the physics of this system. The calculated parameters appear to be compatible with the available data—for example the $`B`$-parameter, relating to the surface energy of the islands shows a similar variation with growth parameters for each island type, which is unsurprising if both island types are dominated by the same facet. However, the fits are in no way perfect, and this may be due to the influence of kinetic effects on the island distributions. Further experiments on the annealing of dot arrays at growth temperature might clarify to what extent the growth proceeds near equilibrium or kinetic effects are predominant.

###### Acknowledgements.

The authors would like to thank O. Briot (Université Montpellier II, France) for providing the GaN templates and H. Mariette (Université Joseph Fourier, Grenoble, France) for many stimulating discussions. R.E.R.’s contribution to this work was performed under the auspices of the U.S. Department of Energy by the University of California, Lawrence Livermore National Laboratory, under Contract No. W-7405-Eng-48.


format: text

text:
# Swapping and entangling qubits of single photons and atoms

## Abstract

A scheme is proposed here to achieve swapping and entangling of photonic and atomic qubits with high fidelity. The mechanism is based on the scattering of a single photon from a $`\mathrm{\Lambda }`$-type three-level atom. The evolution of the coupled system is analyzed by projecting the quantum state onto a ‘bright’ and a ‘dark’ state. Quantum interference of these two states, which is determined by a frequency-dependent phase angle, can be exploited to perform various two-qubit transformations. It is remarkable that the probability of success of such transformations can approach unity in the strong coupling cavity QED regime.

Quantum communication relies on the ability of transmitting and retrieving quantum information at specified locations. Particularly for a quantum network, communication between spatially separated nodes requires inter-conversions between flying qubits and stationary qubits with high fidelity. Therefore it is important to investigate mechanisms that can accomplish the task efficiently. Single photons and long-lived trapped atoms are considered as fundamental hardware in distributed quantum computing, with the former serving ideally as data bus, while the latter playing the role of local quantum memory Monroe-nature . Indeed, various authors have proposed protocols based on strong atom-field coupling in high-$`Q`$ cavities Kimble-transfer ; Haroche-memory . These investigations suggest that a complete transfer of a qubit from an atom to a quantized electromagnetic field can in principle be feasible. Recently, there are also studies addressing methods of quantum communication using ensembles of atoms via their interactions with light in free space Fleischhauer ; Communication-atoms-coherent-light .

In addition to the task of quantum state transfer, it is often desirable to exchange quantum information among carriers of distinct nature Polzik-swapping . A basic operation of this kind is defined by the two-qubit transformation:

$`\left(\alpha _1|1_A+\alpha _0|0_A\right)\left(c_1|1_F+c_0|0_F\right)`$
$`\left(c_1|1_A+c_0|0_A\right)\left(\alpha _1|1_F+\alpha _0|0_F\right),`$ (1)
where $`\alpha _1|1_A+\alpha _0|0_A`$ represents an atomic qubit and $`c_1|1_F+c_0|0_F`$ denotes a photonic qubit. Through the transformation (1), the atomic qubit is deposited to the photon, and at the same time the photonic qubit is mapped onto the atom. Such a swapping process is more general than the typical one-way quantum state transfer previously discussed in cavity QED systems, because the latter corresponds to the special cases with either $`c_1=0`$ (atom to photon) or $`\alpha _1=0`$ (photon to atom) Kimble-transfer ; Haroche-memory . To our knowledge, physical examples of swapping effects have only been investigated in atomic ensembles with collective spin variables Polzik-swapping  and NMR systems nmr , a generic protocol of qubit swapping among single photons and atoms in leaky optical cavities has not yet been found.

In this Letter we present a mechanism to achieve qubit swapping based on scattering of a single photon by an atom in a high-$`Q`$ cavity. The two ground states of a trapped $`\mathrm{\Lambda }`$-type three-level atom play the role of the two states of a stationary qubit, and the two polarizations of a single photon constitute a flying qubit. Quantum interference between the ‘bright’ state and the ‘dark’ state of the $`\mathrm{\Lambda }`$-system holds the key to realization of qubit swapping. We derive explicitly an exact analytic form of the scattering matrix, and determine a frequency-dependent phase shift of the bright state. Such a phase shift controls the quantum interference, and can be large in the strong cavity coupling regime. Therefore the frequency degree of freedom of photon would enable various kinds of two-qubit transformations. In addition to the swapping operation (1), we will also show how a maximally entangled photon-atom pair can be generated at suitable input photon frequencies.

To begin with, we consider a one-dimensional cavity with length $`l`$ and bounded by two mirrors. The left mirror at $`x=0`$ is perfectly reflecting, while the other at $`x=l`$ is partially transparent. The normal modes of the electromagnetic field are characterized by a continuous wave number $`k`$. Specifically, the spatial mode functions $`u_k(x)`$ are given by loudon

$$u_k(x)=\{\begin{array}{cccccccccccccccccccc}I(k)\mathrm{sin}kx& & & & & & & & & & & & & & & & & & & \\ e^{ikx}+R(k)e^{ikx}& & & & & & & & & & & & & & & & & & & \end{array}\begin{array}{cccccccccccccccccccc}& & & & & & & & & & & & & & & & & & & \end{array}\begin{array}{cccccccccccccccccccc}0& <& x& <& l& & & & & & & & & & & & & & & \\ l& <& x& <& \mathrm{}& & & & & & & & & & & & & & & \end{array},$$
(2)
where $`I(k)=2it/(1+re^{2ikl})`$ and $`R(k)=(rt+re^{2ikl})/(1+re^{2ikl})`$, with $`r`$ and $`t`$ being the reflection and transmission coefficients of the right mirror, respectively mirror . These continuous field modes provide a basis for the field quantization. A $`\mathrm{\Lambda }`$-type three-level atom is located near the center of the cavity. The atom has two degenerate ground states $`|L`$, $`|R`$ and an excited state $`|e`$. The ground state $`|L`$ ($`|R`$) can be excited to $`|e`$ by absorbing a $`|k_L`$ ($`|k_R`$) mode photon as shown in Fig. 1, where the subscripts $`L`$, $`R`$ denote the polarization of the photon. In our scheme, there is one and only one photon involved in the scattering process, and no external classical pump fields are required. To facilitate our discussion, we first neglect the spontaneous emission from the atomic excited state into side modes of the cavity. The loss due to spontaneous decay will be discussed later in the paper.

The Hamiltonian of our model (in units of $`\mathrm{}=c=1`$) is given by

$`\widehat{H}`$ $`=`$ $`\omega _e|ee|+{\displaystyle _0^{\mathrm{}}}𝑑k{\displaystyle \underset{\mu =L,R}{}}k\widehat{a}_{k\mu }^{}\widehat{a}_{k\mu }`$ (3)
$`+{\displaystyle _0^{\mathrm{}}}𝑑k{\displaystyle \underset{\mu =L,R}{}}g_\mu (k)\widehat{a}_{k\mu }|e\mu |+\mathrm{h}.\mathrm{c}..`$
Here, $`\widehat{a}_{k\mu }`$ and $`\widehat{a}_{k\mu }^{}`$ are the annihilation and creation operators associated with the field mode $`u_k(x)`$ with the polarization $`\mu `$ ($`\mu =L,R`$). The dipole coupling strength $`g_\mu (k)`$ is given by

$$g_\mu (k)=\frac{\lambda _\mu \sqrt{\kappa /\pi }e^{i\theta _\mu }}{kk_c+i\kappa },$$
(4)
which is proportional to the mode strength at the location of the atom, i.e., $`u_k(x=l/2)`$. The $`\kappa =\mathrm{ln}|r|/2l`$ is the leakage rate of the cavity, $`\lambda _\mu ^2=_{\mathrm{}}^{\mathrm{}}|g_\mu (k)|^2𝑑k`$, and $`\theta _\mu `$ is a phase angle associated with dipole transition matrix elements. In writing Eq. (4), we have assumed that only one of the cavity quasi-modes (with a resonance frequency $`k_c`$ defined by the real part of the pole of $`R(k)`$) interacts with the atom. Such an approximation is valid when $`\omega _ek_c`$, and all other quasi-modes are far off resonance. This requires a small cavity with high finesse in typical cavity QED experiments Kimble-logic-gate . Note also that in the optical regime, the approximation allows us to extend the lower bound of the frequency of the photon from $`0`$ to $`\mathrm{}`$.

Initially, the atom is prepared in a coherent superposition of the two ground states and a single photon wave packet of arbitrary spectrum and polarizations is injected into the cavity. Our task is to determine the analytic solution of the final state of the system. First, we notice that the scattering process is energy conserving. Hence for a given $`k`$, the initial and final states dwell in the same $`k`$-subspace spanned by the four basis vectors: $`|L;k_L`$, $`|L;k_R`$, $`|R;k_L`$, and $`|R;k_R`$. It is obvious that $`|L;k_R`$ and $`|R;k_L`$ are eigenvectors of $`\widehat{H}`$. In addition, the coherent superposition $`g_R(k)|L;k_Lg_L(k)|R;k_R`$ is a dark state well known in $`\mathrm{\Lambda }`$-systems. The excited state $`|e`$ couples only with the ‘bright’ state

$$|\psi (k)\frac{1}{V(k)}\left[g_L^{}(k)|L;k_L+g_R^{}(k)|R;k_R\right].$$
(5)
where $`V(k)=\sqrt{|g_L(k)|^2+|g_R(k)|^2}`$. Therefore we only need to solve the evolution of $`|\psi (k)`$, which is governed by the Hamiltonian (3) in the corresponding subspace:

$`\widehat{H}^{}`$ $`=`$ $`\omega _e|e;\varphi e;\varphi |+{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑kk|\psi (k)\psi (k)|`$ (6)
$`+{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑kV(k)|e;\varphi \psi (k)|+\mathrm{h}.\mathrm{c}.,`$
with $`\varphi `$ being the vacuum of the field.

For the scattering process considered here, we are interested in the transition matrix element $`U_{kk^{}}=\psi (k)|\widehat{U}(t)|\psi (k^{})`$ in the asymptotic long time limit, where $`\widehat{U}(t)`$ is the evolution operator. With the help of the resolvent formalism Cohen , we obtain the relation:

$$U_{kk^{}}=\delta (kk^{})\left[12\pi iV(k)^2e;\varphi |\widehat{G}(k)|e;\varphi \right]e^{ikt},$$
(7)
where $`\widehat{G}(k)(k\widehat{H}^{})^1`$ is the resolvent. The element $`e;\varphi |\widehat{G}(k)|e;\varphi `$ can be determined by the projection operator method, which gives

$$e;\varphi |\widehat{G}(k)|e;\varphi =\frac{\mathrm{\Delta }k+i\kappa }{(\mathrm{\Delta }k\omega _+)(\mathrm{\Delta }k\omega _{})}.$$
(8)
Here $`\mathrm{\Delta }k=kk_c`$, $`\omega _\pm =(\delta _ei\kappa )/2\pm \sqrt{\left(\delta _e+i\kappa \right)^2/4+\lambda _L^2+\lambda _R^2}`$ is related to vacuum Rabi-splitting Kimble\_splitting , and $`\delta _e=\omega _ek_c`$. Together with Eq. (7), we have $`\widehat{U}(t,0)|\psi (k)=e^{ikt}e^{i\delta _s(k)}|\psi (k)`$. The exact expression of the phase shift $`\delta _s(k)`$ is given by

$$e^{i\delta _s(k)}=\frac{(\mathrm{\Delta }k\delta _e)(\mathrm{\Delta }k^2+\kappa ^2)(\mathrm{\Delta }k+i\kappa )(\lambda _L^2+\lambda _R^2)}{(\mathrm{\Delta }k\delta _e)(\mathrm{\Delta }k^2+\kappa ^2)(\mathrm{\Delta }ki\kappa )(\lambda _L^2+\lambda _R^2)}.$$
(9)
Let us now denote a general input state (at a given $`k`$) by $`|\mathrm{in};k=\alpha _1(k)|L;k_L+\alpha _2(k)|R;k_R+\alpha _3(k)|L;k_R+\alpha _4(k)|R;k_L`$, and the corresponding output state by $`|\mathrm{out};k=\beta _1(k)|L;k_L+\beta _2(k)|R;k_R+\beta _3(k)|L;k_R+\beta _4(k)|R;k_L`$. With the help of Eq. (9) and taking account of the (non-interacting) dark state, the input-output relation can be conveniently expressed in a matrix equation,

$$\left(\begin{array}{c}\beta _1\\ \beta _2\\ \beta _3\\ \beta _4\end{array}\right)=\left(\begin{array}{cccc}T_{LL}& T_{LR}& 0& 0\\ T_{RL}& T_{RR}& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right)\left(\begin{array}{c}\alpha _1\\ \alpha _2\\ \alpha _3\\ \alpha _4\end{array}\right).$$
(10)
Here $`T_{LR}(k)=T_{RL}(k)e^{2i(\theta _R\theta _L)}=g_L^{}(k)g_R(k)(e^{i\delta _s(k)}1)/V(k)^2`$, $`T_{LL}(k)=e^{i\delta _s(k)}(|g_L(k)|^2+|g_R(k)|^2)/V(k)^2`$, $`T_{RR}(k)=e^{i\delta _s(k)}(|g_R(k)|^2+|g_L(k)|^2)/V(k)^2`$, and we have omitted the trivial free evolution phase factor $`e^{ikt}`$.

Equation (10) describes a general transformation of photon-atom states before and after the scattering process. The matrix elements $`T_{ij}(k)`$ $`(i,j=L,R)`$ are the consequences of interference between amplitudes associated with the dark state and the bright state $`|\psi (k)`$ that acquires a phase shift $`\delta _s(k)`$. Notice that the elements $`T_{ij}(k)`$ form a sub-matrix that is unitary. By tuning the photon frequency and cavity parameters, we are able to perform qubit transformations between the atom and the photon. In the following, we consider swapping and entangling of qubits.

Swapping qubits— If we treat the labels $`L(k_L)`$ and $`R(k_R)`$ as logical $`1(0)`$ and $`0(1)`$ for the atom(photon), respectively, the qubit swapping operation (1) corresponds to the conditions: $`T_{LL}=T_{RR}=0`$ and $`T_{LR}=T_{RL}=1`$. Such conditions are satisfied if $`g_L(k)=g_R(k)`$ and

$$\delta _s(k)=\pm \pi .$$
(11)
The condition $`g_L(k)=g_R(k)`$ can be satisfied by choosing suitable atomic transition schemes. For example, the D1 line of sodium with hyperfine ground states $`|L=|F=1,m_F=1`$ and $`|R=|F=1,m_F=1`$, and the excited state $`|e=|F=1,m_F=0`$ is a candidate of $`\mathrm{\Lambda }`$-systems with equal but opposite dipole matrix elements.

If the cavity is tuned at the resonance $`k_c=\omega _e`$, the requirement of a $`\pi `$ phase shift $`\delta _s(k)`$ is satisfied for the photon frequencies:

$$k=k_c,k_c\pm \sqrt{2\lambda ^2\kappa ^2},$$
(12)
where $`\lambda \lambda _L=\lambda _R`$. Hence, if the frequency of the incident photon is given by Eq. (12), perfect qubit swapping can be achieved as long as the loss due to spontaneous decay is ignorable.

Furthermore, we can extend our analysis to incorporate a non-zero spontaneous decay rate $`\gamma `$ of the excited state. This can be done by adding a negative imaginary part $`i\gamma `$ to the atomic frequency $`\omega _e`$. Accordingly, Eq. (9) is modified but the general form of input-output transformation given by Eqs. (10) remains unchanged. We find that among the three characteristic frequencies given in Eq. (12), $`k=k_c`$ provides a robust performance against spontaneous emission loss in the strong coupling regime, where $`\lambda ^2\kappa \gamma `$ and consequently $`e^{i\delta _s(k_c)}1+\kappa \gamma /\lambda ^2`$. Therefore the loss due to spontaneous decay is of order $`\kappa \gamma /\lambda ^2`$, and becomes insignificant in this regime.

To complete our analysis we need to consider realistic photons in forms of wavepackets, instead of being purely monochromatic light. We consider an initial state $`|\mathrm{\Phi }_{\mathrm{in}}=\left(A_L|L+A_R|R\right)_{\mathrm{}}^{\mathrm{}}𝑑kf_S(k)\left[C_L|k_L+C_R|k_R\right]`$. To apply our solution (7)-(10), the spectral function $`f_S(k)`$ of the photon packet should be taken such that the atom does not experience the field of the injected photon for $`t0`$. For concreteness, we consider a Gaussian photon packet which is traveling towards the cavity from a far distance $`x_0`$ at $`t=0`$ with a peak frequency $`k_c`$ and spectral width $`\kappa _{\mathrm{in}}`$ ourpra :

$$f_S(k^{})=\frac{1}{\pi ^{1/4}\sqrt{\kappa _{\mathrm{in}}}}\mathrm{exp}\left[\frac{(k^{}k_c)^2}{2\kappa _{\mathrm{in}}^2}+ik^{}x_0\right].$$
(13)
In the asymptotic long time limit, the output state $`|\mathrm{\Phi }_{\mathrm{out}}`$ is determined by the transformation (10) for each $`k`$. The swapping fidelity is defined by the overlap $`F=|\mathrm{\Phi }_{\mathrm{swap}}|\mathrm{\Phi }_{\mathrm{out}}|^2`$, where $`|\mathrm{\Phi }_{\mathrm{swap}}=\left(C_R|L+C_L|R\right)_{\mathrm{}}^{\mathrm{}}𝑑kf_S(k)\left[A_R|k_L+A_L|k_R\right]`$ is the ideal swapping of the input state.

After some straightforward calculations, we have $`F=12\mathrm{R}\mathrm{e}(\xi )\eta +|\xi |^2\eta ^2`$, where $`\eta =|A_LC_LA_RC_R|^2`$ and $`\xi =_{\mathrm{}}^{\mathrm{}}T_{LL}(k)|f_S(k)|^2𝑑k`$. Therefore the fidelity depends on two parameters $`\eta `$ and $`\xi `$, which in turn depend on the qubit states and the photon spectrum. In essence $`\eta `$ measures the overlap of the initial state with the interacting state $`|\psi (k)`$ defined in Eq. (5), and $`\xi `$ measures the deviations of $`T_{ij}`$ from the ideal swapping matrix elements $`T_{LL}=0`$ and $`T_{LR}=1`$ averaged over the incident photon spectrum.

Further calculations show that $`F`$ attains minimum when $`\eta =1`$ for any $`\xi `$, therefore $`F`$ obeys the inequality,

$$F12\mathrm{R}\mathrm{e}(\xi )+|\xi |^2=F_{\mathrm{min}}.$$
(14)
To provide numerical examples, the circles in Fig. 2 show the dependence of $`F_{\mathrm{min}}`$ on the normalized coupling strength $`\lambda /\kappa `$. In the case $`\lambda /\kappa =10`$, $`\gamma /\kappa =0.5`$, we have $`F_{\mathrm{min}}97\%`$ for a Gaussian packet $`\kappa _{in}=0.1\kappa `$. Higher $`F_{\mathrm{min}}`$ can be achieved by decreasing the spectral width of the input photon.

Entangling qubits — If $`T_{RL}(k)=T_{LL}(k)e^{i\theta }`$ (where $`\theta `$ is a real phase angle) is satisfied at certain frequencies, then an initial product state $`|L,k_L`$ at those frequencies will evolve into a Bell’s state: $`|\mathrm{\Phi }_E=(|L,k_L+e^{i\theta }|R,k_R)/\sqrt{2}`$. In other words, the scattered photon and the final state of the atom can become maximally entangled.

We find that solutions of $`k`$ satisfying $`T_{RL}(k)=T_{LL}(k)e^{i\theta }`$ do exist. Under the condition $`\delta _e=0`$, we have $`\mathrm{\Delta }k=\pm \kappa `$ (or $`k=k_c\pm \kappa `$) as approximate solutions in the strong coupling regime $`\lambda \gamma ,\kappa `$, and the corresponding $`\theta `$ are $`\pi /2`$. Such photon frequencies do not depend on the coupling strengths, hence the entangling operation is not sensitive to the position of the atom as long as the strong coupling requirement is satisfied.

To discuss the success probability, let us consider an initial state $`|\mathrm{\Phi }_{\mathrm{in}}=|L_{\mathrm{}}^{\mathrm{}}𝑑kf_E(k)|k_L`$, where $`f_E(k)`$ peaks at $`k=k_c+\kappa `$ with a spectral width $`\kappa _{\mathrm{in}}`$. The success probability is defined by $`P=\left|\mathrm{\Phi }_E|\mathrm{\Phi }_{\mathrm{out}}\right|^2`$, where $`|\mathrm{\Phi }_E=_{\mathrm{}}^{\mathrm{}}𝑑kf_E(k)(|L,k_Li|R,k_R)/\sqrt{2}`$ is the target Bell’s state. It can be shown that: $`P=\frac{1}{2}+|\xi |^2\mathrm{Re}(\xi )+\mathrm{Im}(\xi )`$, where $`\xi `$ is similarly defined by the previous expression, with $`f_S(k)`$ replaced by $`f_E(k)`$. It is useful to remark that $`P`$ depends on the photon spectrum $`|f_E(k)|^2`$. The phase of $`f_E(k)`$, which defines the spatial details of the photon wavefunction, does not affect the success probability. The same is also true for the swapping process above. We find that the probability of generating the Bell’s state is quite high in the strong coupling regime as long as the spectral width of the incident photon is sufficiently narrow. Figure 2 illustrates how $`P`$ depends on the normalized coupling strength $`\lambda /\kappa `$. For example at $`\lambda /\kappa =10`$, $`P99.2\%`$ can be achieved by using $`\gamma =0.5\kappa `$ and a Gaussian amplitudes $`f_E(k)`$ with $`\kappa _{\mathrm{in}}=0.1\kappa `$.

To conclude, we have demonstrated how a single photon and an atom can effectively ‘communicate’ with each other through scattering in a strong coupling cavity QED environment. An interesting aspect of our approach is that the frequency of the photon serves as a control parameter. Various qubit transformations can be achieved in a tunable way according to the scattering matrix (10). The continuous dependence of frequencies in such matrix not only allows swapping and entangling in the same setup, but also enables the preparation of a wide range of photon-atom states when combining local transformations of the atom and the photon individually. Although we have used a one-dimensional cavity to illustrate the process, the Hamiltonian (3) itself is quite general. Our method is equally applicable to systems in higher dimensions as long as a single quasi-mode and a dominant input-output channel is involved in the process. In addition, since our protocol concerns only asymptotic states, no precise timing of interaction period is required as compared with typical state control strategies Swap\_remark . Finally, we note that the realization of our scheme relies on deterministic single photon sources Law-Kimble ; Rempe-photon\_source  and trapping of single atoms inside optical cavities Mundt ; Mckeever . Our work in fact addresses a novel application of single photons scattering once these technologies become available.

###### Acknowledgements.

We thank H.T. Fung for useful discussions. Our work is supported in part by the Hong Kong Research Grants Council (grant No. CUHK4282/00P and CUHK4016/03P).


text:
# Coherence of trapped one-dimensional (quasi-)condensates and continuous atom lasers in waveguides

## I Introduction

An atom laser is a device which emits an intense, directed coherent atomic beam in much the same way as a conventional (optical) laser emits a beam of coherent photons. Interest in such ‘devices’, whose potential applications include fields such as precision metrology, atom holography and nanotechnology, arose after the experimental realization of Bose-Einstein condensation (BEC) in trapped atomic alkali gases . Although BEC need not necessarily be a prerequisite for the production of an atom laser, it was soon realized that a coherent atomic beam could be generated by the coherent extraction of a trapped BEC, thus leading to the first demonstration of a ‘rudimentary atom laser beam’ in 1997 . In this pioneering experiment, atoms from a BEC trapped in a particular magnetic state were coherently extracted by means of pulsed radio-frequency (rf) transitions to a magnetically untrapped state; the outcoupled atoms then fell under gravity. Since then, a number of experimental groups have demonstrated the ability to produce pulsed , or quasi-continuous  atom lasers, with atoms usually extracted from the magnetically trapped state via rf or Raman transitions (although gravity-induced tunnelling from an optically trapped BEC has also been demonstrated as an outcoupling mechanism). The outcoupled atoms propagate away from the centre of the trap with an accelerated motion due to gravity, or at a constant speed due to the momentum imparted to the atoms by the two-photon Raman transitions. Raman outcoupling has the additional advantage of offering good spatial selectivity within the condensate. The removal of trapped atoms depletes the condensate and currently sets a very tight constraint on the duration of atom laser operation. Significant steps to overcome this limitation were recently carried out at MIT, where it was demonstrated that the time required for the creation of a reasonably-sized BEC can be significantly smaller than its lifetime . This has enabled researchers to make a second BEC, merge it with an existing one, thus demonstrating the first real-time replenishing of a condensate. However, at the time of writing this paper, it is still unkown how such merging affects the coherence properties of the original condensate. All above discussion focuses on the production of an atom laser beam based on the coherent extraction a of ground-state BEC, and an interesting topic meriting further discussion is the possibility of creation of coherent atomic beams with different spatial modes, based on the coherent extraction of a non-ground-state BEC . Finally, one should mention an alternative approach currently underway experimentally , which aims to create a continuous atom laser by evaporative cooling in a magnetic waveguide .

The coherence properties of trapped condensates have been studied experimentally using a variety of techniques. Spatial (equilibrium) properties have been studied by interfering two independent BEC’s , by studying two-  and three-body interactions , by spectroscopic techniques , and by interfering atoms extracted from different points within a single condensate , or at different times . The temporal coherence of an atom laser beam extracted from a pre-formed BEC has been measured in , thus also yielding an upper limit for temporal phase fluctuations within the condensate. On the theoretical side, the coherence properties of atom lasers have been studied in numerous publications, and an excellent review of the early literature on this topic (e.g. based on rate equations ) can be found in . Our theoretical understanding of these issues has further improved following the successful analysis of a recent experiment by means of the Gross-Pitaevskii equation  and an in-depth discussion of output coupling based on Raman transitions . However, such treatments are aimed at describing the coherence properties of an atom laser as it is extracted from a pre-formed (pure) condensate. These treatments hence suffer from two limitations: Firstly, they completely ignore the effect of thermal atoms which are inevitably present in the system and can lead to additional phase diffusion of the atom laser beam . Perhaps more significantly for our current discussion, such treatments do not consider a re-pumping mechanism into the condensate, to compensate for the outcoupled atoms. These treatments are therefore not well suited for describing steady-state atom laser operation. The latter necessitates a mechanism pumping atoms into the condensate, and this is usually modelled phenomenologically . In particular, in recent discussions , both trapped and outcoupled components are described by coupled nonlinear Schrödinger equations, with the trapped component additionally pumped from a thermal cloud, in such a manner, that the condensate is continually replenished during the outcoupling.

However, such pumping is assumed to be coherent, such that any atom falling into the magnetic trap becomes immediately part of the condensate, whereas we believe that a realistic pumping mechanism would, in general, introduce an additional phase diffusion of the condensate, as a result of the newly-pumped atoms. To account for this additional feature, one must combine the description of condensate formation based on quantum kinetic theory, with the process of coherent outcoupling. It is the aim of this paper to develop such a description and perform a qualitative and quantitative analysis of the coherence properties of a continuous one-dimensional atom laser. The results presented here are somewhat preliminary and will be further investigated in a subsequent publication. The important new feature of our model is that pumping arises naturally from the non-equilibrium treatment of the coupled dynamics of condensed and thermal atoms in a trap . This gives rise to a Langevin equation for the growth of a condensate in contact with a heat bath  (note that such an approach differs considerably from previous stochastic treatments applied to similar issues ). The outcoupled atoms are then described by a nonlinear Schrödinger equation, which is coupled to the stochastic condensate field by (externally-induced) Raman coupling. The fluctuations inherent in the trapped component are then passed onto the atom laser beam, whose coherence properties we can thus determine.

This work focuses entirely on one dimensional systems, which are known to suffer from enhanced phase fluctuations and possess no off-diagonal long-range order. In particular, if such phase variations occur on length scales considerably smaller that the size of the system, then the system is known as a quasi-condensate . This significant difference in the coherence properties of one-dimensional systems from the three-dimensional case usually encountered in current atom laser experiments indicates that our approach will yield estimates for the atom laser coherence which will be strongly dependent on the coherence properties of the original trapped (quasi-)condensate, and will thus, in general, not be directly applicable to current experiments. For this reason, our current work does not focus on the ‘absolute’ coherence of the atom laser, but rather on the dependence of the coherence properties of the output beam on the outcoupling procedure (given a particular coherence of the trapped system), and the change in the coherence properties of the trapped component due to the outcoupling. To some extent, such conclusions should also hold qualitatively for the three-dimensional systems. In this paper we discuss a continuous raman output coupling scheme to a magnetically untrapped state, with the atoms leaving the ‘interaction region’ at the centre of the trap by means of the momentum imparted to them by the Raman lasers . Throughout this work we assume that the entire system (condensate plus atom laser) remains kinematically one-dimensional, which is equivalent to an atom laser propagating in a transversely very tightly confining waveguide . For technical reasons, the preliminary results presented in this paper are limited to the region of quasi-condensation, with parameters closely resembling those of a recent one-dimensional experiment . The case of (almost) true condensation, and a more in-depth theoretical interpretation of the predictions of our model, will be discussed in a separate publication.

This paper is structured as follows. Sec. 2 discusses in detail the theoretical framework used to study the coherence properties of a (quasi-)condensate in a trap and the atom laser beam. In Sec. 3 we present results for the coherence length of the trapped component, prior to the outcoupling (Sec. 3.1) and after the outcoupling has been turned on (Sec. 3.2). Sec. 4 discusses the coherence length for the outcoupled component, and compares its behaviour at steady state to that of the trapped (quasi-)condensate. In particular, we discuss explicitly how the coherence length of the atom laser beam is modified by changing either the magnitude, or the phase of the external electromagnetic field (corresponding to the momentum imparted to the outcoupled atoms) producing the outcoupling. Finally, Sec. 5 presents some concluding remarks.

## II Theoretical Model

In this section we present the theoretical framework to be used in the rest of the paper. Firstly, we describe the one-dimensional (quasi-)condensate in the trap.

### A Trapped Component Only

Consider an ultracold magnetically trapped three-dimensional gas above the critical point, on which a laser beam is applied in such a manner that it provides an additional tightly-confining optical potential along two of the directions, in a manner similar to a recent experiment . The laser beam can be arranged in such a manner that the potential ‘dimple’ it creates transversely to the waveguide has only one bound state below the chemical potential of the three-dimensional thermal cloud. In this case, the motion of the system becomes ‘frozen out’ along the transverse directions. The system thus becomes kinematically one-dimensional, but remains in contact with the three-dimensional thermal cloud, which acts as its heat bath. The dynamics of the order parameter $`\mathrm{\Phi }(z,t)`$ in the dimple is thus governed by

$$i\mathrm{}\frac{\mathrm{\Phi }(z,t)}{t}=\left[\frac{\mathrm{}^2^2}{2m}+V^{\mathrm{ext}}(z)\mu iR(z,t)+g|\mathrm{\Phi }(z,t)|^2\right]\mathrm{\Phi }(z,t)+\eta (z,t),$$
(1)
where the external trapping potential in the weakly-confining direction is given by $`V^{\mathrm{ext}}(z)=m\omega _z^2z^2/2`$, and $`\mu `$ is the effective chemical potential of the one-dimensional system. The one-dimensional coupling constant $`g`$ is given by $`g=4\pi \mathrm{}^2\kappa /m`$, where $`\kappa `$ corresponds to a one-dimensional scattering length; this is related to the usual three-dimensional scattering length $`a_{3D}`$ by $`\kappa =a_{3D}/2\pi l_{}^2`$, where $`l_{}=\sqrt{(\mathrm{}/m\omega _{})}`$ is the harmonic oscillator length of the axially symmetric trap in the direction perpendicular to the waveguide axis z. Note that $`\mathrm{\Phi }(z,t)`$ corresponds to the entire field in the dimple, and not simply to the condensed component. The role of the thermal component is evident in the contributions $`iR(z,t)`$ and $`\eta (z,t)`$. Physically, the function $`iR(z,t)`$ describes the pumping of the one-dimensional gas from the surrounding reservoir, and $`\eta (z,t)`$ corresponds to the associated noise with Gaussian correlations. In the classical field approximation  (which is well-satisfied in current experiments), these are given by

$$iR(z,t)=\frac{\beta }{4}\mathrm{}\mathrm{\Sigma }^\mathrm{K}(z)\left(\frac{\mathrm{}^2^2}{2m}+V^{\mathrm{ext}}(z)\mu +g|\mathrm{\Phi }(z,t)|^2\right),$$
(2)
$$\eta ^{}(z,t)\eta (z^{},t^{})=\frac{i\mathrm{}^2}{2}\mathrm{\Sigma }^\mathrm{K}(z)\delta (zz^{})\delta (tt^{}),$$
(3)
where $`\mathrm{}`$ denotes averaging over the realizations of the noise $`\eta (z,t)`$, and $`\beta =1/(k_BT)`$ as usual. Both above quantities depend on the Keldysh self-energy $`\mathrm{}\mathrm{\Sigma }^K(z)`$ which arises physically from collisions that scatter an atom out of, or into the ‘heat bath’, and is given by

$`\mathrm{}\mathrm{\Sigma }^K(z)`$ $`=`$ $`4i{\displaystyle \frac{g^2}{(2\pi )^5\mathrm{}^6}}{\displaystyle 𝑑p_1𝑑p_2𝑑p_3\delta (p_1p_2p_3)\delta (ϵ_1ϵ_2ϵ_3+V_{ext})}`$ (5)
$`\times (N(ϵ_1)+1)N(ϵ_2)N(ϵ_3)`$
where $`ϵ_i=(p_i^2/2m)+V_{ext}(z_i)`$ and $`N(ϵ_i)=[exp(\beta (ϵ_i\mu ))1]^1`$ is the Bose occupation factor. The numerical techniques employed are discussed in Ref. , where it was also shown that the above expressions guarantee that the trapped gas relaxes to the correct equilibrium, as ensured by the fluctuation-dissipation theorem. To simplify the numerics, the noncondensed part in the dimple is here allowed to relax to the “classical” value $`N(ϵ)=[\beta (ϵ\mu )]^1`$.

Fig. 1 shows the evolution of the total density profile in the one-dimensional trap at different relaxation times (after turning on the optical dimple), for a sodium (quasi-)condensate at $`T=200nK`$ and with trapping frequencies taken from a recent experiment . The effect of the thermal contributions is obvious, since the equilibrium density profile differs considerably from the corresponding (zero-temperature) Thomas-Fermi solution obtained in the absence of noise.

### B Coupled Evolution of Trapped and Outcoupled Components

Having discussed how a one-dimensional (quasi-)condensate is produced, we next describe our model for a one-dimensional atom laser. After allowing the trapped component to relax to its equilibrium value, the outcoupling mechanism is turned on, leading to a continuous depletion of the (quasi-)condensate. Nonetheless, the presence of the additional pumping mechanism from the surrounding thermal cloud (which acts as its heat bath), leads to a steady-state population of the trapped component. An implicit crucial assumption made here is that the heat bath remains unaffected as it pumps atoms into the (quasi-)condensate, to compensate for the atoms lost from the trap due to outcoupling. This would indeed be the case (at least for short times) in the limit of very slow outcoupling and large initial thermal cloud. In our model, we envisage that the ultracold thermal cloud is continually replenished at a rate comparable to that of outcoupling. Then, our model is expected to be valid as long as the residual oscillations in the heat bath atom number (which indirectly affect the coherence of the trapped component and of the atom laser) are kept arbitrarily small, which we believe lies within current experimental capabilities. As noted earlier, we further assume that our entire system remains kinematically one-dimensional, which requires the optical dimple to extend over the entire range of our atom laser.

To model the raman outcoupling scheme, we follow the usual procedure of coupling the order parameter of the trapped system to a nonlinear Schrödinger equation for the outcoupled component. Such an approach has been already discussed in the literature , by making the further crucial assumption that the trapped component can be adequately described by a nonlinear Schrödinger equation with phenomenological pumping. Our treatment improves on such earlier attempts in that we explicitly include the effect of the thermal cloud on the traped component by means of a non-equilibrium theory based on the Schwinger-Keldysh formalism . Hence, in our treatment pumping arises as a result of the relaxation of the system into its new confining potential.

The coupled evolution of the trapped $`\mathrm{\Phi }_1(z,t)`$ and outcoupled $`\mathrm{\Phi }_2(z,t)`$ components is thus determined by

$`i\mathrm{}{\displaystyle \frac{\mathrm{\Phi }_1(z,t)}{t}}`$ $`=`$ $`\left[{\displaystyle \frac{\mathrm{}^2^2}{2m}}+V_1^{\mathrm{ext}}(z)\mu +\delta iR(z,t)+g_{11}|\mathrm{\Phi }_1(z,t)|^2+g_{12}|\mathrm{\Phi }_2(z,t)|^2\right]\mathrm{\Phi }_1(z,t)`$ (7)
$`+\mathrm{\Omega }(z,t)\mathrm{\Phi }_2(z,t)+\eta (z,t)`$
$`i\mathrm{}{\displaystyle \frac{\mathrm{\Phi }_2(z,t)}{t}}`$ $`=`$ $`\left[{\displaystyle \frac{\mathrm{}^2^2}{2m}}+V_2^{\mathrm{ext}}(z)\mu +g_{22}|\mathrm{\Phi }_2(z,t)|^2+g_{12}|\mathrm{\Phi }_1(z,t)|^2\right]\mathrm{\Phi }_2(z,t)`$ (9)
$`+\mathrm{\Omega }^{}(z,t)\mathrm{\Phi }_1(z,t).`$
Here $`g_{ii}`$ correspond to the one-dimensional self-interactions of each component, whereas $`g_{12}`$ is the mean-field interaction between the two components. For simplicity, we have taken $`g_{11}=g_{22}=g_{12}=g`$ (defined earlier). The two components are coupled by the application of an external electromagnetic field via $`\mathrm{\Omega }(z,t)=\mathrm{\Omega }_0e^{ikz}`$, where $`\mathrm{\Omega }_0`$ the strength of the Raman coupling and $`\mathrm{}k`$ the momentum imparted to the outcoupled atoms. Our description is limited to a one-dimensional atom laser in a horizontal waveguide, and thus ignores gravity. For $`{}_{}{}^{23}Na`$, the trapped atoms are assumed to be in the $`|F,m_F=|1,1`$ state (for which $`V_1^{\mathrm{ext}}(z)=m\omega _z^2z^2/2`$), whereas for simplicity we have chosen to outcouple the atoms to the untrapped $`|F,m_F=|1,0`$ state (for which $`V_2^{\mathrm{ext}}(z)=0`$). All our simulations are performed at zero detuning ($`\delta =0`$) such that the outcoupling, which is peaked precisely at $`z=0`$, occurs only in a very narrow region around the centre of the trap. An illustrative discussion of the effect of changing the detuning can be found in .

Fig. 2 shows the density profiles of the trapped and outcoupled components at steady state, for two different coupling parameters $`\mathrm{\Omega }_0`$. We see clearly that the momentum kick imparted to the outcoupled atoms, leads to their propagation in the negative z-direction (note that absorbing boundaries have been used at the edges of our computational grid). Due to the mean field interaction between trapped and outcoupled components, such directional motion of the outcoupled atoms creates an asymmetry in the profile of the trapped atoms, which is more pronounced in the case of large $`\mathrm{\Omega }_0`$. Since the outcoupling is primarily from the trap centre ($`\delta =0`$), this leads to a “double-peak” structure for the trapped atom density profile.

Having discussed our model and its predictions for the density profiles we now turn to the discussion of the coherence properties, which we investigate in terms of an appropriately defined coherence length.

### C Definition of Coherence Length

The coherence properties of the system can be inferred from the behaviour of the normalized first-order off-diagonal correlation function defined by

$$g^{(1)}(z_0,t_0;z,t)=\frac{\mathrm{\Phi }_\alpha ^{}(z_0,t_0)\mathrm{\Phi }_\alpha (z,t)}{\sqrt{|\mathrm{\Phi }_\alpha (z_0,t_0)|^2|\mathrm{\Phi }_\alpha (z,t)|^2}},$$
(10)
where $`\mathrm{}`$ denotes averaging over noise realizations, and $`\alpha =1,2`$ denote respectively the trapped and outcoupled components. Such a quantity can be calculated by numerical autocorrelation measurements. Fig. 3 shows the position dependence of the equal-time correlation function $`g^{(1)}(z_0,t_0;z,t_0)`$ for the trapped component at the trap centre ($`z_0=0`$), once it has relaxed to its equilibrium value in the absence of outcoupling. The observed rapid decay of the correlation function over the size of the system is clear evidence of quasi-condensation, as discussed in . Since this closely resembles an exponential , we have chosen to define the coherence length, $`l_{coh}`$ as the value of $`z`$ for which this function decays to $`(1/e)`$ (more generally $`g^{(1)}(z_0,t_0;(z_0+l_{coh}),t_0)=(1/e)`$. In this paper we will investigate the behaviour of $`l_{coh}`$ (i) as a function of $`z_0`$ for both trapped and outcoupled components, at fixed relaxation time $`t_0`$, and (ii) as a function of relaxation time $`t_0`$, for fixed position $`z_0`$. Note that by relaxation time we define the time elapsed from the moment that the tightly confining transverse potential is turned on.

## III Coherence Length of Trapped System

### A Without Outcoupling

Initially we discuss the relaxation of the trapped component to its equilibrium value in the dimple prior to the outcoupling. Fig. 4(a) shows the temporal evolution of the coherence length at various positions in the trap. It is seen clearly that $`l_{coh}`$ grows significantly at the trap centre and remains mostly unaffected at the edges of the trap (where there are only thermal atoms present). The dependence of $`l_{coh}`$ on position in the trap at various times can also be seen in Fig. 4(b), which confirms the anticipated result that the coherence length is largest at the trap centre (where the quasi-condensate density is maximum).

### B With Outcoupling

The coherence length of the trapped component can decrease dramatically once the outcoupling is turned on, as a result of the removal of coherent atoms from the trap centre where the Raman transition resonance condition is satisfied. The values of the coherence length of the trapped component once steady state is reached between pumping and outcoupling are (potentially much) lower than the original values prior to outcoupling, because the atoms pumped from the surrounding thermal cloud to compensate for the atoms lost due to the outcoupling mechanism, do not necessarily share the same phase as the trapped atoms before their outcoupling. This is shown clearly in Fig. 5(a), which also demonstrates that the decrease in coherence is more dramatic at the trap centre (since atoms are mostly removed from that region) and essentially negligible outside the Thomas-Fermi radius of the system. Fig. 5(a) focuses on the regime of strong outcoupling where these effects are pronounced. In a typical experiment, one might seek to minimize the perturbations on the trapped component due to the outcoupling, and such a case (corresponding to the ‘weak outcoupling’ of Fig. 2(b)) will be further considered in Sec. 4.

Fig. 5(b) shows the dependence of the coherence length of the trapped component on position, once steady state has been reached. Apart from the rapid decrease of the coherence length due to the outcoupling procedure, we find that the coherence length is no longer maximum at the trap centre, but this maximum in now slightly offset along the z-axis. This reflects the fact that the maximum of the quasi-condensate density is slightly offset from its original value (as seen in Fig. 2; note that the positive z-axis of Figs. 3-9 corresponds to the negative z-axis of Fig. 2), due to the asymmetry created from the mean-field interaction between trapped and outcoupled atoms as a result of the directional motion of the outcoupled component.

## IV Coherence Length of Atom Laser Beam

In order to discuss the coherence of an atom laser, we now perform the same type of analysis for the outcoupled component. Fig. 6 shows the coherence length of the atom laser beam as a function of time for particular points along the beam, chosen such that they lie well outside the Thomas-Fermi region where the interaction with the trapped quasi-condensate significantly affects its behaviour. From this figure, we deduce that the atom laser is essentially operating under steady-state conditions. Fig. 7 shows the position dependence of the coherence length of the trapped and outcoupled components at steady state, over the entire range of the trap. The coherence length of the outcoupled component is found to change considerably within the region where the trapped component exhibits significant coherence (ie essentially within the Thomas-Fermi radius of the quasi-condensate), and reaches a steady-state value far away from the center of the trap. On the contrary, the coherence length of the trapped component (apart from the minor increase at small values of $`|z|`$ due to the asymmetry in the quasi-condensate density profile discussed above) decreases significantly as we move away from the centre of the trap. This shows that the outcoupled component acts as an atom laser, in the sense that it has a reasonably constant coherence length both in space and in time.

All figures shown so far focused on the case of strong outcoupling which heavily depletes the quasi-condensate (since, in this limit, effects become more noticeable). On the contrary, the situation of most experimental interest is likely to be that in which the trapped component is only slightly perturbed by the outcoupling. This regime is discussed in Fig. 8 which plots the position dependence of the coherence length of both quasi-condensate and atom laser, for two different coupling strengths $`\mathrm{\Omega }_0`$, focusing on the region of interest sufficiently far from the trap centre. The coherence length of the trapped component decreases rapidly towards the edges of the trap (since there is no off-diagonal long range order in one dimension) and is larger in the case of weaker outcoupling, since this leads to less depletion of the quasi-condensate (with the effect more pronounced within the Thomas-Fermi region which is not shown here). We believe that it is for the same reasons, that we find a longer steady state atom laser coherence length in the case of stronger outcoupling (see also ).

Finally, Fig. 9 investigates the effect of imparting different ‘momentum kicks’ $`k`$ to the outcoupled atoms, at a constant outcoupling rate $`\mathrm{\Omega }_0`$. The coherence length of the trapped component in the region $`z_0<12l_z`$ is found to be somewhat larger in the limit of smaller momentum kicks. A possible interpretation, may be the following: since the outcoupled atoms leave the central ‘interaction region’ at a slower rate, the non-markovian nature of the mean field interactions between the two components  ensures that there is always a larger fraction of trapped atoms in this case, and hence more coherence within the trapped component. On the other hand, the coherence length of the atom laser beam is significantly increased for larger momenta, due to the increased speed of propagation of the outcoupled atoms in the waveguide.

## V Conclusions

The coherence properties of an atom laser beam are still only partially understood, despite significant work done on this field both experimentally and theoretically over the past few years. In this paper, we have focused on a one-dimensional atom laser beam in a tightly-confining waveguide, and we have shown that it can indeed maintain a reasonably constant coherence length in both space and time, in the region sufficiently far away from the centre of the trap (where the perturbations due to interactions with the trapped coherent component are minimized). Our system was assumed to be kinematically one-dimensional, and thus our analysis completely suppressed the effect of the transverse directions of the waveguide in which the atom laser propagates. In three-dimensional systems, the effects of the transverse modes can become significant, as discussed in . Our values for the coherence length of the atom laser beam are considerably smaller than the corresponding three-dimensional systems, because it is well-known that phase fluctuations in low-dimensional systems can be strongly enhanced, leading to the appearance of quasi-condensation (that is, condensation with rapidly varying phase over the size of the system), as opposed to true condensation which arises in three-dimensional systems. In principle, this ‘defect’ can be cured, by working in that region of the phase diagram, where one essentially obtains ‘true BEC’ . A further limitation of the approach discussed in this paper is that for computational reasons, the thermal cloud was treated clasically, which means that its size was slightly overestimated, leading to a minor underestimate in the coherence properties of our system. Nonetheless, simulations of the full quantum theory have shown this effect to be a rather minor correction, at least for the parameters considered here .

Notwithstanding the above remarks, we consider our current contribution as yet another important step in the description of the coherence properties of atom lasers since, to the best of our knowledge, our work is the first qualitative and quantitative comparison of the coherence length of (quasi-)condensates and atom lasers beyond the usual mean field treatments. The analysis presented in this paper is based on somewhat preliminary results, and we intend to further investigate and support our conclusions in a subsequent publication. Furthermore, we are currently investigating the related issues of coherence time and linewidth of such atom lasers for experimentally relevant parameters (and in the limit of almost pure condensation), and such results will be presented elsewhere.

ACKNOWLEDGEMENTS

I am indebted to Professor Henk Stoof for suggesting the application of the particular stochastic field method to study the coherence properties of atom lasers, and for his invaluable comments during the early stages of this work. I would also like to acknowledge computational assistance from Mark Leadbeater, and the hospitality of the Benasque Center for Science, where I had the opportunity to discuss my results with various participants of the ‘Physics of Ultracold Dilute Atomic Gases’ Workshop. Funding for this work was provided by the UK EPSRC.


text:
# Large-scale electronic structure calculation and its application

## I Introduction

Nanoscale materials are directly governed by quantum mechanical freedoms of electron systems. Nowadays, electron systems of realistic materials are treated by the ab initio electronic structure calculations that are based on the density functional theory (DFT) HOHENBERG-KOHN ; KOHN-SHAM , the first-principle molecular dynamics CP , and related theories developed for decades. A typical system size of present ab initio calculations is, however, on the order of 102 atoms and new practical theories are required for nanoscale calculations. This article is devoted to the methods in large-scale electronic structure calculations and their application to nanoscale materials Hoshi00a ; Hoshi01 ; Hoshi02b ; GESHI ; TAKAYAMA ; Hoshi03a .

In general, a quantum mechanical calculation of an electron system is reduced to an eigen value equation;

$`\widehat{H}|\varphi _k^{(\mathrm{eig})}=\epsilon _k^{(\mathrm{eig})}|\varphi _k^{(\mathrm{eig})}`$ (1)
with an effective one-body Hamiltonian $`\widehat{H}`$. Here the eigen energies and eigen states are denoted as $`\{\epsilon _k^{(\mathrm{eig})}\}`$ and $`\{\varphi _k^{(\mathrm{eig})}\}`$, respectively. A physical quantity $`\widehat{X}`$ is given as

$`\widehat{X}={\displaystyle \underset{k}{\overset{\mathrm{occ}.}{}}}\varphi _k^{(\mathrm{eig})}|\widehat{X}|\varphi _k^{(\mathrm{eig})}=\mathrm{Tr}[\widehat{\rho }\widehat{X}]`$ (2)
with occupied eigen states $`\{\varphi _k^{(\mathrm{eig})}\}`$ or the one-body density matrix $`\widehat{\rho }`$

$`\widehat{\rho }{\displaystyle \underset{k}{\overset{\mathrm{occ}.}{}}}|\varphi _k^{(\mathrm{eig})}\varphi _k^{(\mathrm{eig})}|.`$ (3)
The calculation of eigen states, Eq. (1), is usually reduced to a matrix diagonalization procedure and gives a severe computational cost. Therefore, the essential methodology for large-scale calculations is how to obtain the density matrix $`\widehat{\rho }`$ without calculating eigen states. This article focuses the methods for structural properties including molecular dynamics simulations. For the above purpose, the most import physical quantity is the total energy NOTE .

The trace in Eq. (2) is expressed as

$`\mathrm{Tr}[\widehat{\rho }\widehat{X}]={\displaystyle 𝑑𝒓𝑑𝒓^{}\rho (𝒓,𝒓^{})X(𝒓^{},𝒓)}`$ (4)
with real space coordinates $`𝒓`$,$`𝒓^{}`$. The off-diagonal components of the density matrix, $`\rho (𝒓,𝒓^{}),𝒓𝒓^{}`$, are essential quantum mechanical freedoms, while the diagonal components, the charge density at the point $`𝒓`$ $`(\rho (𝒓,𝒓)n(𝒓))`$, appear also in classical mechanics. As an important fact for practical large-scale calculations, the off-diagonal long range component of the density matrix dose not contribute explicitly to the value of $`\widehat{X}`$, if the operator $`\widehat{X}`$ is a short range one. We can found a general principle within DFT, called ‘nearsightedness principle’KOHN96 , which is directly related to the above fact.

For practical algorithms, there are many proposals. See reviews or comparison papers ORDEJON-REV ; Goedecker99a ; GALLI-REV ; JAYANTHI-REV ; BOWLER97A . In this article, we pick out two methods; (i) method with generalized Wannier state and (ii) Krylov subspace method. Figure 1 demonstrates the computational cost with diagonalization and our calculationHoshi02b . Here one can find that the diagonalization results in a compuational cost proportional to $`N^3`$ with the system size $`(N)`$, as is usual in matrix diagonalization procedure ($`N^3`$). Our calculation, on the other hand, shows an ‘order-$`N`$’ property, with upto 106 atoms, in the sense that the computational cost is proportional to the system size ($`N`$).

## II Generalized Wannier state

The generalized Wannier state is a generalization of the (conventional) Wannier states WANNIER-ORG ; KOHN-WANI59 ; NOTE-AM  (See Appendix A). Its pioneering works were done by Walter Kohn in the context of large-scale calculations KOHN-WANI73 ; KOHN-WANI93 . The pictures of the generalized Wannier states are localized ‘chemical’ wave functions in condensed matter, such as a bonding orbital or a lone-pair orbital, with a slight spatial extension or ‘tail’.

The generalized Wannier states $`\{\varphi _i^{(\mathrm{WS})}\}`$ are defined as localized wave functions that satisfy the equation

$`H|\varphi _i^{(\mathrm{WS})}={\displaystyle \underset{j=1}{\overset{\mathrm{occ}}{}}}\epsilon _{ij}|\varphi _j^{(\mathrm{WS})}`$ (5)
and the orthogonality

$`\varphi _i^{(\mathrm{WS})}|\varphi _j^{(\mathrm{WS})}=\delta _{ij}.`$ (6)
The matrix $`\epsilon _{ij}`$ is introduced as the Lagrange multiplier for the constraint of Eq. (6) and is given as

$`\epsilon _{ij}=\varphi _j^{(\mathrm{WS})}|H|\varphi _i^{(\mathrm{WS})}.`$ (7)
The solutions of Eq. (5) is equivalent to the unitary transformation of the eigen states $`\{\varphi _k^{(\mathrm{eig})}\}`$

$`|\varphi _i^{(\mathrm{WS})}={\displaystyle \underset{k}{\overset{\mathrm{occ}.}{}}}U_{ik}|\varphi _k^{(\mathrm{eig})},`$ (8)
where $`U_{ik}`$ is a unitary matrix. Here the suffix $`i`$ of the Wannier state $`\varphi _i^{(\mathrm{WS})}`$ denotes its localization center.

It is crucial that the generalized Wannier states reproduce the one-body density matrix $`\widehat{\rho }`$ in Eq. (3), where the eigen states $`\{\varphi _k^{(\mathrm{eig})}\}`$ are replaced by the Wannier states $`\{\varphi _j^{(\mathrm{WS})}\}`$. In results, any physical quantity can be reproduced in the trace form of Eq. (2). NOTE-WS-LMTO .

The concept of the generalized Wannier state is used for practical large-scale calculations LOM-MGC ; LOM-MGC2 ; Hoshi00a . We derived a mapped eigen value equation for the generalized Wannier states

$`H_{\mathrm{WS}}^{(i)}|\varphi _i^{(\mathrm{WS})}=\epsilon _{\mathrm{WS}}^{(i)}|\varphi _i^{(\mathrm{WS})}`$ (9)
with a mapped Hamiltonian $`H_{\mathrm{WS}}^{(i)}`$ that is dependent on the other Wannier states $`\{\varphi _j^{(\mathrm{WS})}\}_{ji}`$ Hoshi00a ; Hoshi01 . Equation (9) is equivalent to Eqs. (5) and (6). A Wannier state $`|\varphi _i^{(\mathrm{WS})}`$ is not an eigen state of the original Hamiltonian $`H`$ but an eigen state of the above mapped Hamiltonian $`H_{\mathrm{WS}}^{(i)}`$. Equation (9) also shows that the locality of a Wannier state can be mapped, formally, to that of a virtual impurity state Hoshi00a .

With Eq. (9), we developed a variational method so as to generate approximate Wannier states, which is called variational Wannier state method Hoshi00a . As the practical procedure, Eq. (9) is solved iteratively under explicit localization constraint on each Wannier state

$`\{\varphi _i^{(\mathrm{WS})}\}\{H_{\mathrm{WS}}^{(i)}\}\{\varphi _i^{(\mathrm{WS})}\}\{H_{\mathrm{WS}}^{(i)}\}\mathrm{}`$ (10)
With Eq. (9), we also developed a perturbative method to generate Wannier states, which is called perturbative Wannier state method Hoshi00a ; Hoshi01 . This method corresponds to a non-iterative solution of Eq. (9). It is noteworthy that the perturbative Wannier state, unlike the variational one, is localized without any explicit localization constraint, when a short range Hamiltonian $`H`$ is used.

## III Krylov subspace method

The Krylov subspace method KRYLOV  is focused as another important concept for large-scale calculations. The Krylov subspace is a mathematical concept and its definition is the linear space that is constructed from the following vectors;

$$|i,H|i,H^2|i,\mathrm{}H^{\nu 1}|i.$$
(11)
Here an initial vector $`|i`$ should be given. In the present context, the matrix $`H`$ is a Hamiltonian. The number of bases in the Krylov subspace ($`\nu `$) is chosen to be much smaller than that of the original Hamiltonian matrix $`H`$. In a practical method for large-scale calculations, we consider the Hamiltonian matrix only within the above subspace, which means the drastic reduction of the matrix size. The Krylov subspace gives the mathematical foundation of many numerical algorithms such as the standard conjugate gradient methodKRYLOV  and the recursion methodHAYDOCK-REV ; PETTIFOR89 .

Recently, we developed a practical Krylov subspace method for the calculation of the density matrix and applied it to molecular dynamics simulations TAKAYAMA . The number of bases in the subspace was chosen, typically, as $`\nu =30`$. We compared the resultant density matrix with that by the Wannier state method. Now we are planning other molecular dynamics simulations, especially, to metals.

## IV Hybrid scheme and parallel computations

As another fundamental methodology for large-scale calculations, we developed the hybrid scheme within quantum mechanics Hoshi02b ; Hoshi03a . The basic idea is the followings; The one-body density matrix is decomposed into two partial matrices or ‘subsystems’ that are constructed from several occupied wave functions. This decomposition corresponds to dividing the occupied Hilbert space. The different partial density matrices are solved by different solver methods. Each subsystem is obtained with a well-defined mapped Hamiltonian and a well-defined electron numberHoshi03a . Test calculations are done by the combinations between (a) the diagonalization method and perturbative Wannier state methods, (b) the variational and perturbative Wannier state methods Hoshi02b  (c) the Krylov subspace method and the perturbative Wannier state method Hoshi03a . Since the present hybrid scheme is a technique in calculating the density matrix $`\widehat{\rho }`$, any physical quantity is quantum mechanically well defined with Eq. (2).

Parallel computation is also important for large-scale calculations. Test calculation of the perturbative Wannier state method is carried out with upto 106 atoms GESHI  using the Message Passing Interface techniqueMPI-URL  and with upto 107 atoms Hoshi03a  using the OpenMP technique OMP-URL . We are now developing the parallelization of other methods TAKAYAMA ; Hoshi03a .

## V Application and discussion

As a practical nanoscale application, the molecular dynamics simulation is performed for fracture of nanocrystalline silicon Hoshi02b . A standard workstation is used for the simulations with upto 105 atoms. We use the hybrid scheme between the variational Wannier state method and the perturbative Wannier state method.

In the continuum theory of fracture GRIFFITH ; NOTE2 , a critical crack length is defined by a dimensional analysis between the competitive energy terms of the bulk strain (3D) energy and the surface formation (2D) energy. Since the definition of the critical length is independent on the sample size or the lattice constant, we can expect a crossover in fracture phenomena between nanoscale and macroscale samples. The investigation of the above crossover is one purpose of the present simulation. Another purpose is the fracture behavior in atomistic pictures, on the points of how and why the fracture path is formed and propagates in the crystalline geometry NOTE-FRAC .

Dynamical fracture processes are simulated under external loads in the  direction. As the elementary process in fracture, we observe a two-stage surface reconstruction process. The process contains the drastic change of the Wannier states from the bulk (sp3) bonding state to surface ones. Figure 2 shows a result, in which the fracture propagates anisotropically on the (001) plane and reconstructed surfaces appear with asymmetric dimers Hoshi02b . Step structures are formed in larger systems so as to reduce the anisotropic surface strain energy within a flat (001) surface. Such a step formation is understood as the beginning of a crossover between nanoscale and macroscale samples Hoshi02b . Further investigation should be done for direct discussion of the crossover.

The present calculations are carried out using tight-binding Hamiltonian within $`s`$ and $`p`$ orbitals. We should say that its applicability is rather limited, due to the simplicity of Hamiltonian. Its parameter theory, however, reproduces systematically several ab initio results among different elements or phases, because the tight-binding formulation is universal within the scaled length and energy units Hoshi00a ; Hoshi01 ; Hoshi03a . An important future work is to construct simple and practical (tight-binding) Hamiltonians more systematically from the ab initio theory. We will use the muffin-tin orbital formulation for the construction, because it gives directly the tight-binding formulation LMTO-HIGHLIGHT .

Recently the concept ‘multiscale mechanics’ is focused as the seamless theoretical connection of material simulation methods among the three principles of mechanics; (I) quantum mechanics (for electron systems), (II) classical mechanics, and (III) continuum mechanics. The present work gives a guiding principle and a typical example for the concept, which is carried out by simplifying the total energy functional.

## Appendix A Derivation of Conventional Wannier state

Here we derive the conventional Wannier state WANNIER-ORG ; NOTE-AM  as a specific case of Eq. (8). In periodic systems, eigen states are called Bloch states $`\{\psi _{\nu 𝒌}^{(\mathrm{Bloch})}\}`$ with the suffices of the band $`\nu `$ and the k-point $`𝒌`$, the point in the Brillouin zone. Within an isolated single band, the Wannier states can be defined $`W_{\nu 𝒍}`$ with the suffices of the band $`\nu `$ and the lattice vector $`𝒍`$;

$`W_{\nu 𝒍}(𝒓)={\displaystyle 𝑑𝒌e^{i𝒌𝒓}\psi _{\nu 𝒌}^{(\mathrm{Bloch})}(𝒓)},`$ (12)
where the integration is done within the Brillouin zone NOTE4 . Equation (8) will be reduced to Eq. (12) , when the corresponding unitary matrix $`U`$ is chosen as

$`U_{ij}U_{\nu 𝒍,\nu ^{}𝒌}\delta _{\nu \nu ^{}}e^{i𝒌𝒓}.`$ (13)
The conventional Wannier state is given by the unitary transform only within an isolated single band $`(\nu =\nu ^{})`$, while the generalized Wannier states are given by the unitary transform within different bands $`(\nu \nu ^{})`$. It is also noteworthy that the concept of the generalized Wannier state, unlike the conventional one, can be applicable to non-periodic cases.