text:
# Angle-resolved photoemission study and first principles calculation of the electronic structure of GaTe

## I Introduction

Layered materials of the III-VI family like GaSe, InSe, and GaTe are of special interest for their potential applications in nonlinear and optical bistable devices, as well as in the development of solar cells and solid-state batteries. Among these compounds, GaTe is one of the less studied because it presents a higher anisotropic and more complex crystal structure than the rest of the III-VI materials. GaTe belongs to the B2/m space group and has a base centered monoclinic unit cell. The layer structure \[see Fig. 1(a)\], as that of GaSe or InSe, is composed of a four-sheet anion-cation-cation-anion intralayer stacking pattern in which the bonds are mainly covalent with some ionic contribution. Nevertheless, in contrast to GaSe and InSe, one-third of the cation-cation bonds in GaTe lies almost in the layer plane in a direction perpendicular to the c axis. Differences also exist when considering layer stacking sequences. Similarly to the rest of III-VI compounds, interlayer anion-anion bonds in GaTe are of weak van der Waals type. Generally, this fact promotes the existence of polytypes with different stacking sequences, as is observed in InSe and GaSe. Nevertheless, no polytypism seems to occur in GaTe, as is suggested by the absence of conjugate modes in Raman-scattering and infrared-absorption measurements. The whole structure gives rise to fourfold coordination for the Ga atoms, with three Te and one Ga atoms, and to threefold coordination for the Te with Ga atoms. As a consequence of these differences in the crystal structure of GaTe with respect to other layered compounds, GaTe has only one twofold symmetry axis, lying along the c axis, whereas other III-VI compounds have a threefold symmetry axis perpendicular to the layer plane.

Besides the common in-plane-out-of-plane anisotropy observed in the optical and electrical properties of the layered compounds, the special in-plane anisotropy in bonding in GaTe makes the optical properties also anisotropic in the layer plane. GaTe is a direct-gap semiconductor with strong excitonic absorption even at room temperature. The optical properties near the band-gap energy region appear to be slightly anisotropic. Anisotropy of the absorption coefficient has been observed to increase at higher energies, which was attributed to transitions from deep states related to the in-plane Ga-Ga bonds. Moreover, low-temperature photoluminescence and transmission measurements revealed the presence of different exciton states depending on the in-plane polarization vector. Recently, these results were explained by the existence of three types of exciton states due to j-j coupling, whose optical transitions are allowed under the group symmetry of this material.

There is still a lack of both theoretical and experimental studies on the electronic band structure of GaTe. One of the reasons for this lacuna lies on the complexity of the primitive unit cell, that contains 6 GaTe units (54 valence electrons), which has made band structure calculations hardly available up to now. Recently, the electronic band structure of GaTe has been calculated along some directions by the ab initio tight-binding linear muffin-tin orbitals method. These results are in agreement with the direct-gap behavior of this semiconductor and locate the band edge at the border of the Brillouin zone (BZ). Nevertheless, there are still some open questions regarding the electronic band structure along the high-symmetry directions of the GaTe BZ, as well as the experimental confirmation of the electronic band structure obtained in the calculations.

In this work we study the electronic properties of the GaTe compound. The GaTe band structure has been calculated by numerical atomic orbitals density-functional theory (NAO-DFT). In addition, angle-resolved photoemission (ARPES) measurements have been carried out in GaTe(001) single crystals. From these measurements, the valence-band dispersion of GaTe was extracted and interpreted on the basis of the band structure calculations. The experimental setup is described in Sec. II. Section III is devoted to the band structure of GaTe calculated along different directions of the GaTe BZ, including the technical aspects of the calculation. In Sec. IV, the experimental results obtained by ARPES are presented and discussed in the frame of the calculated GaTe band structure. Finally, in Sec. V we present the main conclusions of this work.

## II Experimental

The experiments were performed at LURE (Orsay, France) using the Spanish-French (PES2) experimental station of the Super-Aco storage ring, described elsewhere. The measurements were carried out in a purpose-built ultra-high vacuum system, with a base pressure of 5x10-11 mbar, equipped with an angle resolving 50 mm hemispherical VSW analyzer coupled on a goniometer inside the chamber. The manipulator was mounted in a two-axes goniometer which allows rotation of the sample in the whole 360o azimuthal angle and in the 180o polar emission angle relative to surface normal ($`\mathrm{\Theta }_{off}`$), with an overall angular resolution of 0.5o. Photoelectrons were excited with p-polarized synchrotron radiation in the 20-60 eV energy range. For this selected photon energy (h$`\nu `$) range, the energy resolution was of 60 meV. With this setup, different energy-distribution-curve (EDC) series were recorded scanning the h$`\nu `$ as well as $`\mathrm{\Theta }_{off}`$, at constant incident angle of the light ($`\mathrm{\Theta }_i`$). In these conditions, changes of polarization effects on initial states along one series are neglected. All the EDC spectra shown in this work are referred to the Fermi level (EF).

The GaTe single crystals used in this work were cut from an ingot grown by the Bridgman-Stockbarger method. The samples were cleaved in situ after introduction in the ultra-high vacuum chamber. The samples were easily cleaved in the layer plane due to the existence of weak interlayer van der Waals bonds. LEED spectra showed a sharp spot pattern corresponding to the bulk monoclinic material. No surface impurities were detected by synchrotron radiation photoemission measurements. The samples were oriented by azimuthal and polar photoelectron diffraction scans recording the Ga 3d peak intensity. Unambiguous azimuthal orientation of the samples was possible due to the B2/m symmetry of the crystal.

## III Band structure of GaTe

### A Calculation details

In this work, we present fully self-consistent density-functional theory (DFT) calculations of the electronic structure of GaTe. The calculations are performed in the local-density approximation. The exchange-correlation potential is that of Ceperley and Adler as parameterized by Perdew and Zunger. Only the valence electrons are considered in the calculation, with the core being replaced by norm-conserving scalar relativistic pseudopotentials factorized in the Kleinman-Bylander form. The pseudopotentials were generated using the atomic valence configurations 5s25p4 for Te and 3d104s24p1 for Ga. The cutoff radii were 2.0, 2.0, 3.0, and 3.0 a.u. for the s, p, d, and f components in Te, respectively, and 2.1, 2.3, 1.1 and 2.0 a.u. for the s, p, d, and f components of Ga, respectively. Since the core and valence charges overlap significantly for both atomic species, we include nonlinear partial-core corrections with matching radii of 1.0 a.u. for Te and 0.5 a.u. for Ga to describe the exchange and correlation in the core region.

The valence one-particle problem was solved using a linear combination of numerical (pseudo-) atomic orbitals with finite range. The details of the basis generation (including multiple-$`\zeta `$ and polarization functions) can be found elsewhere. In this work we have used a split-valence double-$`\zeta `$ basis set with a single shell of polarization orbitals (that is, containing two s shells, two p shells, and one d shell both for Te and for Ga), as obtained with an energy shift of 250 meV and a split norm of 15%. The integrals of the self-consistent terms of the Kohn-Sham Hamiltonian are obtained with the help of a regular real space grid on which the electron density is projected. The Hartree potential is calculated by means of fast Fourier transforms in that grid. The grid spacing is determined by the maximum kinetic energy of the plane waves that can be represented in that grid. In the present work, we used a cutoff of around 300 Ry (which changes slightly with the volume of the unit cell). A regular grid of 5x5x5 k points was used to sample the BZ. We have checked that the results are well converged with respect to the real space grid, the BZ sampling and the range of the numerical atomic orbitals. The calculations were performed using the Siesta code.

In our band structure calculations, we used a centered monoclinic cell with primitive vectors $`\text{a}\text{}=\text{b}+(\text{c}\text{a})/2`$, $`\text{b}\text{}=\text{b}(\text{a}+\text{c})/2`$, and $`\text{c}\text{}=\text{b}`$, where a, b, and c are the repeat vectors of the monoclinic cell \[see Figure 1(a)\]. Note that, with this election, a and b vectors are in the layer plane.

### B Brillouin zone and electronic band structure

Figure 1(a) shows the crystal structure of GaTe, described previously. In this figure, we have indicated the crystallographic axes of the GaTe monoclinic unit cell. The GaTe crystal structure is quite similar to those of GaSe and InSe, except for the existence of cation-cation bonds lying almost in the layer plane, with the same coordination and bonding behavior. The comparison can be also extended to the reciprocal space. Figure 1(b) shows the calculated GaTe BZ, in which the main points are indicated. The relation between the different symmetry directions of the GaTe BZ and the reciprocal lattice vectors of the centered monoclinic unit cell is shown in Table I.

The GaTe BZ can be easily described in terms of the crystal structure. The binary axis corresponds to the $`\mathrm{\Gamma }`$Y high-symmetry direction and a crystal vector perpendicular to the layer plane defines a direction parallel to the $`\mathrm{\Gamma }`$Z high-symmetry direction, in the reciprocal space. The GaTe BZ resembles that of InSe, but contracted along the $`\mathrm{\Gamma }`$X direction. The contraction is due to the presence of the quasi-in-plane Ga-Ga bonds in GaTe, which produces that the fundamental translation in the layer plane consists of three Ga2Te2 unities, in contrast to one unity in the other III-VI compounds. On the other side, as both GaTe and InSe have only one layer per unit cell, one should expect the electronic band structure of GaTe to be related to that of InSe, after a three-folding in the layer plane.

Figure 2 shows the electronic band structure of GaTe along different directions of the GaTe BZ \[Fig. 1(b)\] calculated by the NAO-DFT method described previously. The origin of binding energy has been taken at the EF. In spite of the large number of bands obtained, a detailed description of the valence band can be done. The GaTe unit cell consists of three Ga2Te2 unities. Therefore, there are per unit cell: three Ga-Ga bonds, 18 Ga-Te bonds, and one electronic pair associated to each one of the six Te. This gives rise to 27 doubly occupied bands, which are precisely those appearing with binding energies between 14 eV and EF and composing the GaTe valence band. Below this energy (at 16.5 eV), the 3d-core levels of Ga are present (not shown). The unoccupied side of the band structure is rather more complicated to be described, since, besides the 18 Ga-Te and 3 Ga-Ga antibonding levels, all contributions from the polarization d-orbitals overlap.

A detailed band to band description of the band structure of GaTe would not be very useful due to the fact that the low symmetry of the crystal makes the orbital mixture of each band to be very high. Nevertheless, a general description of the orbital character of the GaTe band structure can be performed by means of the density of states (DOS). Figure 3(a) shows the total DOS calculated by the NAO-DFT method. Figure 3(b) and 3(c) show the different projections of the s, p, and d orbitals onto the DOS for Ga and Te, respectively. These results indicate that the valence band of GaTe can be divided into three different groups: (i) One group of six bands appearing between 14 and 11 eV, which show a pronounced Te 5s relative behavior with a small Ga 4s and 4p contribution, (ii) a second group of six bands appearing between 8 and 5 eV, which mostly show a Ga 4s character with a pronounced Te 5p contribution in the low binding energy side, and (iii) a third group of 15 bands comprising the range of energy between 5 eV and the EF, which are a mixture of Te 5p and Ga 4p orbitals with a small Ga 4d and 4s contribution. The origin of these groups differs from one to the other. On one side, the first and second group of bands appear to be mostly originated by the Ga-Te bonds, showing those with a marked p-like behavior a more pronounced dispersion along the $`\mathrm{\Gamma }`$Z high-symmetry direction (as well as along directions perpendicular to it) than those with a predominant s-like behavior (Fig. 2). On the other side, in the third group of bands, they coexist six bands originated also from the Ga-Te bonds, three bands from the Ga-Ga bonds, and the six electronic pairs associated to the Te (mostly Te 5p). Focusing our attention on the topmost valence bands, the fact that their orbital character has an important p-like component perpendicular to the layers, makes that most of the bands display a noticeable dispersion along $`\mathrm{\Gamma }`$Z high-symmetry direction (Fig. 2).

With regard to the DOS corresponding to the unoccupied side of the band structure, the huge overlap of Ga-Te and Ga-Ga antibonding levels with the d-orbitals of polarization makes an orbital description quite more complicated than that carried out for the occupied side of the DOS. Nevertheless, it should be mentioned that states at the conduction band edge are mostly originated from Te 5p and Ga 4p orbitals, with an important contribution from Ga 4s and a smaller contribution Te 5d and Ga 4d, which suggests a pronounced spz-like behavior to the conduction band edge. Only for higher states in the conduction band, contributions from Ga and Te d-orbitals of polarization appear to be of the same magnitude than those of the p-orbitals.

The general shape and the different orbital contributions to the calculated DOS by NAO-DFT are quite similar to those obtained previously by the ab initio tight binding method. Some differences can be appreciated in the unoccupied side of the DOS. By both methods, the behavior of the conduction band edge seems to be established as a mixture of s and p orbitals. Nevertheless, in contrast to previous calculations, the Ga 4p contribution to the conduction band edge appears to be of the same order as that of the Ga 4s, which gives a more pronounced p-like behavior to the conduction band edge in the present calculations. Also the contribution of Ga and Te d-orbitals of polarization to the conduction band edge, as well as to deeper states of the conduction band, seem to be underestimated in those calculations (mostly the Ga 4d). As regards the band dispersion, authors seem to have used a standard six-face monoclinic BZ, which does not correspond to the GaTe BZ. The band dispersion calculated along some directions in Ref.  (see Table I) is very similar to the present results. However, along other directions, as the $`\mathrm{\Gamma }`$Z and AM directions in Fig. 8 of Ref. , the path crosses out of the first BZ to neighboring ones, giving rise to the quasi-symmetrical band dispersion observed along these directions.

From the above results, it can be concluded that GaTe is a direct-gap semiconductor with the gap located at zone border at Z point of the GaTe BZ, which has the same symmetry as $`\mathrm{\Gamma }`$ (B2/m). The band gap appears to be underestimated (1.098 eV), which is a well-known tendency of the LDA approximation in semiconductors. As in InSe, the valence-band maximum is mostly formed by antibonding pz-orbitals of the anion. Nevertheless, in contrast to InSe, the conduction-band minimum is not only formed by the s-orbital contribution of the cation. In GaTe a more intense contribution of cation and anion p-orbitals can be observed. This can be attributed to the fact that the presence of the quasi-in-plane cation-cation bonds distorts the anion disposition and, therefore, the interaction of Te atoms with those of contiguous layers and with Ga atoms increases.

## IV Experimental Results and discussion

In this section, we compare the calculated electronic band structure of GaTe with the band dispersion determined experimentally by ARPES in GaTe(001) single crystals. First, we will show and discuss the electronic band dispersion along the $`\mathrm{\Gamma }`$Z high-symmetry direction, proceeding later with those determined along different directions perpendicular to it.

Figure 4(a) shows an EDC series recorded at normal emission with constant $`\mathrm{\Theta }_i`$ in a GaTe(001) sample, as a function of the h$`\nu `$. In these spectra, several peaks have been identified and labelled by small solid bars. In some of them, two sharp dispersing peaks can be also observed, which can be attributed to second-order transitions from deeper occupied states induced by the experimental system. As the h$`\nu `$ increases, these sharp peaks move to lower binding energies. Therefore, photoemission signal from deeper valence-band transitions than those of Fig. 4(a) can be studied in EDCs recorded by extending up the h$`\nu `$-range of Fig. 4(a). Figure. 4(b) shows the results of these measurements. In these spectra, only broad features can be identified, at binding energies of 12.6 eV, which have been labelled as s-band.

By these results, the dispersion of states along the $`\mathrm{\Gamma }`$Z high-symmetry direction can be determined. Nevertheless, to compare these results with the calculated band structure, the momentum of photoelectrons should be estimated. Photoemission process involves energy and crystal momentum conservation. In the frame of the three-step model of the photoexcitation mechanism in ARPES, both momentum components, k and k, of the photoelectrons inside a bulk material can be expressed as

$$k_{}=\sqrt{\frac{2m}{\mathrm{}^2}}\sqrt{h\nu BE\mathrm{\Phi }}\mathrm{sin}\left(\mathrm{\Theta }_{off}\right)$$
(2)
and

$$k_{}=\sqrt{\frac{2m}{\mathrm{}^2}}\sqrt{\left(h\nu BE\mathrm{\Phi }\right)\mathrm{cos}^2\left(\mathrm{\Theta }_{off}\right)V_o},$$
(3)
where BE is the binding energy, $`\mathrm{\Phi }`$ is the work function, Vo is the inner potential, m is the free-electron mass, and $`\mathrm{}`$ is the reduced Planck constant. A complete determination of the photoelectron momentum requires an accurate estimate of the Vo \[Eq.(3)\]. We have assumed it as a constant parameter. An initial testing value of Vo$``$-19 eV ($`\mathrm{\Phi }`$$``$6 eV) can be obtained from the bottom of the valence band (Fig. 2). After testing several values of Vo and taking into account that kΓZ=0.421 Å-1, we have adopted a value of Vo=-25.7 eV, which is two times higher than that found for InSe. Such a high value of the Vo for GaTe is not surprising, since the orbital hybridization giving rise to the GaTe valence-band structure is higher than in InSe. In any case, it is reasonably close to the value expected from the valence-band calculations.

Figure 4(c) shows the valence-band dispersion extracted from the peaks identified in Fig. 4(a), by means of Eqs. (IV). The calculated valence-band dispersion along the $`\mathrm{\Gamma }`$Z high-symmetry direction has been also included. In spite of the large number of different bands present, it can be observed that the calculated band structure reproduces most of the traces observed in the experimental band dispersion. The large concentration of quite flat bands appearing at lower binding energies is well reproduced by the calculated valence-band structure, with states lying at similar binding energies. These facts also put forward the highly hybridized behavior of top states of the valence-band structure predicted by calculations (Fig. 2). Besides this, dispersion of states with a marked Te 5p-Ga 4s orbital character (that is, states between 8 and 5 eV in binding energy) is quite well reproduced. Moreover, the presence of the broad s-band peaks in Fig. 4(b) can be attributed to transitions from the Te 5s-like states, located at binding energies of $``$12 eV (Fig. 2). Nevertheless, some discrepancies appear for states lying at the top of the valence band. Experimentally, the valence-band maximum along the $`\mathrm{\Gamma }`$Z high-symmetry direction appears to be at Z point and the minimum of the valence-band lies at $`\mathrm{\Gamma }`$. These facts are in concordance with the calculated valence-band structure. Moreover, dispersion of the topmost band around Z is well reproduced by the calculations. However, the valence-band minimum at $`\mathrm{\Gamma }`$ obtained by calculations appears shifted down by 0.8 eV with respect to the experimental one. A similar situation occurs with the second topmost valence band. For this band, the calculated band minimum, also at $`\mathrm{\Gamma }`$, appears to be shifted down by 0.6 eV with respect to the experimental one. This disagreement will be discussed later on.

Figure 5(a) shows an EDC series measured by ARPES with h$`\nu `$=33 eV in a GaTe(001) sample. The sample orientation has been chosen to probe, by these ARPES measurements, points of the reciprocal space contained in the $`\mathrm{\Gamma }`$ZH plane. In fact, for this selected h$`\nu `$, these points are expected to have a k lying quite close to the Z point \[Ecs.(IV)\]. In the spectra shown in Fig. 5(a), several peaks have been identified and labelled by small solid bars. Figure 5(b) shows the band dispersion extracted from the peaks identified in Fig. 5(a).

In order to compare the experimental band dispersion shown in Fig. 5(b) with the calculated band structure, it should be taken into account that, by the ARPES measurements summarized in Fig. 5(a), particular points of the reciprocal space along a parabolic curve are probed for each constant initial state \[Eqs. (IV)\]. For initial states with low binding energy, the minimum of the parabola will be close to the Z point, but when running from Z to H (kZH=0.744 Å-1) the parabolic curve perpendicularly deviates from horizontality by 0.07 Å-1. For deeper initial states, the parabolas probe points somewhere along the $`\mathrm{\Gamma }`$Z high-symmetry direction. These facts produce that, by these measurements, states along the ZH direction are not strictly probed. Nevertheless, it has been found that dispersion of valence-band states along the $`\mathrm{\Gamma }`$Z high-symmetry direction can be considered negligible in a large k-range \[Fig. 4(c)\]. These facts encourages us to assume that the calculated band structure along the ZH direction is a good approximation for the experimental band dispersion, as shown in Fig. 5(b). In this figure, it has been also introduced the calculated band dispersion of states with binding energies higher than 3 eV along the $`\mathrm{\Gamma }`$Y high-symmetry direction, which is not far from that calculated along the ZH direction and gives a better agreement with the experimental band dispersion of deeper bands.

Figure 6(a) and 6(b) show two different EDC series measured by ARPES with h$`\nu `$=33 eV in a GaTe(001) sample. In these ARPES measurements, a sample orientation has been chosen to probe points of the reciprocal space contained in the $`\mathrm{\Gamma }`$ZM plane. In these spectra, several peaks have been identified and labelled by small solid bars. Figure 6(c) shows the $`\mathrm{\Gamma }`$ZM transversal cut of the GaTe BZ, in the extended zone scheme. The parabolic curve defined by ARPES measurements with h$`\nu `$=33 eV for constant initial states at the EF is also included in this figure. This scheme illustrates the fact that, by both sets of measurements shown in Figs. 6(a) and 6(b), points of the reciprocal space are probed along the MZM’ and M$`\mathrm{\Gamma }`$M’ directions, respectively. Figure 6(d) shows the band dispersion extracted from the peaks identified in Figs. 6(a) and 6(b). Assuming the same hypotheses as those discussed previously, we have also included the calculated band dispersion along the MZM’ and M$`\mathrm{\Gamma }`$M’ directions, which reproduce most of the traces of the band dispersion obtained experimentally. Focusing again attention on the two topmost valence-bands, their dispersion around Z along the MZM’ direction appears to be rather flat, coherently with the behavior expected by NAO-DFT calculations. Nevertheless, dispersion of these bands along the M$`\mathrm{\Gamma }`$M’ direction shows the same behavior as that observed along the $`\mathrm{\Gamma }`$Z high-symmetry direction \[Figs. 4(c)\], with the same separation at $`\mathrm{\Gamma }`$ between the calculated and experimental minimums of these bands.

The above results obtained by ARPES measurements along different directions of the GaTe BZ support the results obtained by NAO-DFT band structure calculations (Fig. 2). Among these results, it should be emphasized the fact that experimental and calculated band structure place the valence-band maximum in GaTe at zone border at Z, similarly as occurs in the InSe compound. Nevertheless, despite the good agreement obtained through the whole reciprocal space studied, there still remains one point of discrepancy related to the two topmost valence-bands at $`\mathrm{\Gamma }`$ \[Figs. 4(c) and 6(d)\], since the binding energy of the minimum of these bands appear to be overestimated by calculations. These facts can be explained since spin-orbit interaction was not considered in the present calculations. In the related compound InSe, the band giving rise to the valence band maximum at Z mostly has a marked pz-orbital character. States of this band have a total orbital moment j=1/2. Without considering spin-orbit interaction, this band disperses down along the $`\mathrm{\Gamma }`$Z direction, showing the minimum at $`\mathrm{\Gamma }`$. Below this band, the next ones are quite flat bands with a px-py character. Spin-orbit interaction splits off px-py states with different j, rising those with j=3/2 by 0.3-0.4 eV, with respect to the valence band maximum. This produces that the minimum of the topmost valence band at $`\mathrm{\Gamma }`$ is determined by the shifted up px-py bands. A more intense spin-orbit coupling can be expected in GaTe, since the nature of its topmost valence bands is the same as those of InSe and its anion has a higher mass. In fact, the spin-orbit induced shift on the valence-band minimum at $`\mathrm{\Gamma }`$ in GaTe may be evaluated. Taking into account that the atomic spin-orbit p3/2-p1/2 splitting is of 589 and 247 meV for Te and Se, respectively, in GaTe one should expect a shift up of the j=3/2 px-py bands at $`\mathrm{\Gamma }`$ of the order of 0.7-0.9 eV, which is close to that observed in Figs. 4(c) and 6(d).

## V Conclusions

The electronic band structure of GaTe single crystal has been calculated by a NAO-DFT method. These results show that, on one side, the deepest valence bands are mostly formed by Ga-Te orbitals with a marked s-like behavior and, on the other side, bands with a low binding energy have a pronounced p-orbital character originated by cation-anion bonds, as well as by anion electronic pairs and cation-cation bonds. Also, the conduction-band structure was determined. It appears to be composed by bands with mainly a p-orbital character, originated from the antibonding Ga-Te and Ga-Ga bonds, intermixed with the d-orbitals of polarization.

The electronic properties of GaTe single crystal have been studied also by ARPES measurements. By these measurements, dispersion of the valence bands along various directions of the GaTe BZ has been determined. The experimental band dispersion obtained appears to be quite well reproduced by the calculated band dispersion along these directions.

These results establish that GaTe is a direct-gap semiconductor with a gap located at zone border, at Z point. The valence-band maximum shows a marked p-like behavior, with a pronounced anion contribution. The conduction band minimum arises from states with a comparable s\- p-cation and p-anion orbital contribution. Spin-orbit interaction appears to specially alter the dispersion and binding energy of states of the topmost valence bands lying at $`\mathrm{\Gamma }`$. Spin-orbit favors the hybridization of the topmost pz-valence band with deeper and flatter px-py bands. Also, it rises the valence-band minimum at $`\mathrm{\Gamma }`$ towards the EF, which appears to be determined by the shifted up j=3/2 px-py bands.

###### Acknowledgements.

This work was financed by the Large Scale Facilities program of the EU to LURE. The authors gratefully acknowledge Dr. J. Avila and Dr. M.C. Asensio for technical and scientific support at LURE. J.F.S.-R. and G.T. acknowledge financial support from the Ministerio de Educación y Cultura and Ministerio de Ciencia y Tecnología of Spain, respectively. Part of this work was supported by DGI-Spain (Project BFM2000-1312-C02-01), Generalitat de Catalunya (Project 1999SGR207) and Fundación Ramón Areces. The computations were carried out using the resources of CESCA and CEPBA, coordinated by C4.


format: text

text:
# Vortex precession in a rotating nonaxisymmetric trapped Bose-Einstein condensate

## I introduction

Since the original experimental observation of Bose-Einstein condensation in dilute atomic gases , the creation and detection of vortices has attracted great interest (for a recent review article, see, for example Ref. ). In particular, one experiment has measured the precession frequency $`\omega _p`$ of a vortex line in a nonrotating stationary spherical condensate  and compared $`\omega _p`$ to various theoretical estimates . Although these theoretical approaches largely predict the same results, perhaps the most transparent and physical picture  relies on a variational Lagrangian method  that focuses directly on the position of the vortex core. For an axisymmetric condensate in rotational equilibrium at a small angular velocity $`\mathrm{\Omega }\omega _{}`$, the original precession frequency $`\omega _p`$ is altered to $`\omega _p(\mathrm{\Omega })=\omega _p\mathrm{\Omega }`$ by the external rotation. The present work describes the nontrivial generalization of the variational Lagrangian method to a large rotating nonaxisymmetric condensate in the Thomas-Fermi (TF) limit. For extreme asymmetry, we may note that the precession frequency $`\omega _p(\mathrm{\Omega })`$ becomes nearly independent of $`\mathrm{\Omega }`$ (for $`\mathrm{\Omega }\omega _{}`$) because the irrotational flow induced by the rotating asymmetry screens the effect of the external rotation.

## II precessional dynamics of a straight vortex

A Bose condensate is described by a macroscopic order parameter (the condensate wave function) $`\mathrm{\Psi }`$, with particle density $`n=|\mathrm{\Psi }|^2`$. Consider a condensate in a nonaxisymmetric harmonic trap potential $`V_{\mathrm{tr}}(𝐫)=\frac{1}{2}M(\omega _x^2x^2+\omega _y^2y^2+\omega _z^2z^2)`$. The condensate experiences both the trap potential $`V_{\mathrm{tr}}`$ and the self-consistent Hartree interaction $`V_H=g|\mathrm{\Psi }|^2=gn`$ arising from the interaction with all the other particles, where $`g=4\pi \mathrm{}^2a/M`$ is the effective interparticle interaction strength and $`a>0`$ is the $`s`$-wave scattering length . At zero temperature, the order parameter obeys the time-dependent Gross-Pitaevskii (GP) equation

$$i\mathrm{}\frac{\mathrm{\Psi }}{t}=\frac{\mathrm{}^2^2}{2M}\mathrm{\Psi }+V_{\mathrm{tr}}\mathrm{\Psi }+g|\mathrm{\Psi }|^2\mathrm{\Psi }.$$
(1)
When the trap rotates with angular velocity $`\mathrm{\Omega }`$ about the $`z`$ axis, the corresponding GP equation in a co-rotating frame acquires an additional term $`\mathrm{\Omega }L_z\mathrm{\Psi }=i\mathrm{}\mathrm{\Omega }𝐫\times \mathbf{}\mathrm{\Psi }`$ on the right-hand side, arising from the transformation to the rotating frame .

Instead of working directly with the time-dependent GP equation (1), it is often preferable to introduce a Lagrangian formalism based on the Lagrangian functional

$$[\mathrm{\Psi }]=𝒯[\mathrm{\Psi }][\mathrm{\Psi }],$$
(2)
where

$$𝒯[\mathrm{\Psi }]=𝑑V\frac{i\mathrm{}}{2}\left(\mathrm{\Psi }^{}\frac{\mathrm{\Psi }}{t}\mathrm{\Psi }\frac{\mathrm{\Psi }^{}}{t}\right)$$
(3)
is the time-dependent part of the Lagrangian that is analogous to the kinetic energy in classical particle mechanics. Similarly, the free energy $`[\mathrm{\Psi }]`$ here plays the role of the potential energy. In a nonrotating system, $``$ is simply the energy functional $``$, which includes the gradient (bending) energy, the trap energy, and the self-consistent Hartree energy. When the system rotates, however, the free energy has an extra term $`=\mathrm{\Omega }L_z`$ , and the appropriate free-energy functional of the system is

$$[\mathrm{\Psi }]=𝑑V\left(\frac{\mathrm{}^2}{2M}|\mathrm{\Psi }|^2+V_{\mathrm{tr}}|\mathrm{\Psi }|^2+\frac{g}{2}|\mathrm{\Psi }|^4+i\mathrm{}\mathrm{\Omega }\mathrm{\Psi }^{}\frac{\mathrm{\Psi }}{\varphi }\right).$$
(4)
Given the Lagrangian, the associated action is the time integral $`𝒮[\mathrm{\Psi }]=_{t_1}^{t_2}𝑑t[\mathrm{\Psi }]`$, and it is straightforward to verify that the action is stationary with respect to small variations of $`\mathrm{\Psi }`$ and $`\mathrm{\Psi }^{}`$ when $`\mathrm{\Psi }`$ obeys the time-dependent GP equation (1). Thus the time-dependent GP equation is the Euler-Lagrange equation for this problem.

This Lagrangian formalism is exact, but it also provides the basis for a powerful approximate variational method. If the wave function $`\mathrm{\Psi }`$ depends on various parameters, the integral in Eq. (3) will involve the first time derivatives of these parameters. Consequently, the resulting $``$ will serve as the effective Lagrangian that provides the corresponding dynamical equations for these parameters. This method has been used to study the quadrupole oscillations of a vortex-free condensate . Here, we use it to study the dynamics of a vortex in a rotating nonaxisymmetric trap.

Specifically, we assume that a straight singly quantized vortex is displaced from the center of the trap with transverse coordinates $`x_0`$ and $`y_0`$; these variables will serve as time-dependent variational parameters that obey Lagrange’s equations . In the TF limit, the vortex induces negligible change in the condensate density, and, in the presence of the rotation, $`|\mathrm{\Psi }|^2`$ is given by

$$g|\mathrm{\Psi }(𝐫)|^2=\stackrel{~}{\mu }\stackrel{~}{V}_{\mathrm{tr}}(𝐫)=\stackrel{~}{\mu }\frac{1}{2}M\left(\stackrel{~}{\omega }_x^2x^2+\stackrel{~}{\omega }_y^2y^2+\omega _z^2z^2\right)$$
(5)
where $`\stackrel{~}{\omega }_x`$ and $`\stackrel{~}{\omega }_y`$ are effective oscillator frequencies that determine the shape of the condensate in the rotating asymmetric trap. Specifically, they are given in terms of a parameter $`\alpha `$:

$$\stackrel{~}{\omega }_x^2=\omega _x^2+\alpha ^22\alpha \mathrm{\Omega },\stackrel{~}{\omega }_y^2=\omega _y^2+\alpha ^2+2\alpha \mathrm{\Omega },$$
(6)
where $`\alpha `$ satisfies a cubic equation

$$2\alpha ^3+\alpha (\omega _x^2+\omega _y^24\mathrm{\Omega }^2)+\mathrm{\Omega }(\omega _x^2\omega _y^2)=0.$$
(7)
The parameter $`\alpha `$ has a simple physical interpretation , for it determines the irrotational flow $`𝐯_0`$ induced by the rotating nonaxisymmetric trap through the velocity potential $`\mathrm{\Phi }_0=\alpha xy`$, with $`𝐯_0=\mathbf{}\mathrm{\Phi }_0`$. As noted below, $`\alpha =\mathrm{\Omega }\left(\stackrel{~}{\omega }_x^2\stackrel{~}{\omega }_y^2\right)/\left(\stackrel{~}{\omega }_x^2+\stackrel{~}{\omega }_y^2\right)`$.

Equation (5) shows that the condensate density has the familiar TF parabolic form

$$|\mathrm{\Psi }|^2=n_0\left(1\frac{x^2}{\stackrel{~}{R}_x^2}\frac{y^2}{\stackrel{~}{R}_y^2}\frac{z^2}{\stackrel{~}{R}_z^2}\right),$$
(8)
if the right-hand side is positive and zero otherwise. Here $`n_0=\stackrel{~}{\mu }/g=M\stackrel{~}{\mu }/4\pi a\mathrm{}^2`$ is the central density, $`\stackrel{~}{R}_i=\sqrt{2\stackrel{~}{\mu }/M\stackrel{~}{\omega }_i^2}`$ (for $`i=x,y`$), $`\stackrel{~}{R}_z=\sqrt{2\stackrel{~}{\mu }/M\omega _z^2}`$ are the TF radii, and $`\stackrel{~}{\mu }=\frac{1}{2}\mathrm{}\stackrel{~}{\omega }_0\left(15Na/\stackrel{~}{d}_0\right)^{2/5}`$ with $`\stackrel{~}{\omega }_0^3=\stackrel{~}{\omega }_x\stackrel{~}{\omega }_y\omega _z`$ and $`\stackrel{~}{d}_0^2=\mathrm{}/M\stackrel{~}{\omega }_0`$. For a small angular velocity $`\mathrm{\Omega }\omega _{}`$ where $`\omega _{}^2\frac{1}{2}(\omega _x^2+\omega _y^2)`$, the renormalization of the oscillator frequencies and TF radii is of order $`\mathrm{\Omega }^2`$ and thus small.

In our variational approach, we assume that the dominant vortex-induced contribution to the Lagrangian arises from the condensate’s superfluid motion and use the following trial wave function

$$\mathrm{\Psi }=|\mathrm{\Psi }|e^{i(S_0+S_1)},$$
(9)
where $`|\mathrm{\Psi }|`$ follows from the TF density profile in Eq. (8). The first term of the phase is taken as $`S_0=(M\alpha /\mathrm{})xy`$, which represents the irrotational motion induced by the rotating nonaxisymmetric trap (note that $`\alpha `$ vanishes if $`\omega _x=\omega _y`$). The second term is taken to be

$$S_1=\mathrm{arctan}\left(\frac{yy_0}{xx_0}\right),$$
(10)
which characterizes the circulating flow around the vortex line. In contrast to the analysis of Refs. , we do not include an image vortex because the form of the TF condensate density ensures that the particle current vanishes automatically at the TF surface.

Variation of the functional (4) yields an Euler-Lagrange equation for $`S_0`$, whose solution gives

$$S_0=\frac{M\mathrm{\Omega }}{\mathrm{}}\left(\frac{\stackrel{~}{\omega }_x^2\stackrel{~}{\omega }_y^2}{\stackrel{~}{\omega }_x^2+\stackrel{~}{\omega }_y^2}\right)xy=\frac{M\mathrm{\Omega }}{\mathrm{}}\left(\frac{\stackrel{~}{R}_x^2\stackrel{~}{R}_y^2}{\stackrel{~}{R}_x^2+\stackrel{~}{R}_y^2}\right)xy.$$
(11)
Note that the time-dependent parameters $`x_0`$ and $`y_0`$ appear only in the phase $`S_1`$, so that

$`_t\mathrm{\Psi }`$ $`=`$ $`i\mathrm{\Psi }_tS_1`$ (12)
$`=`$ $`i\mathrm{\Psi }\dot{𝐫}_0\mathbf{}S_1`$ (13)
$`=`$ $`i\mathrm{\Psi }{\displaystyle \frac{\dot{x}_0(yy_0)+\dot{y}_0(xx_0)}{(xx_0)^2+(yy_0)^2}}.`$ (14)
## III Two-dimensional condensate

For simplicity, we first consider a two-dimensional condensate that is unbounded in the $`z`$ direction (taking $`\omega _z0`$), in which case the various terms in the Lagrangian are interpreted per unit length. To evaluate the time-dependent part (3), it is convenient to introduce dimensionless coordinates scaled with the renormalized TF radii and shift the origin of coordinates to the position of the vortex. In this way, we find

$`𝒯[\mathrm{\Psi }]`$ $`=`$ $`\mathrm{}n_0\stackrel{~}{R}_x^2\stackrel{~}{R}_y^2{\displaystyle d^2r\frac{\dot{x}_0y+\dot{y}_0x}{\stackrel{~}{R}_x^2x^2+\stackrel{~}{R}_y^2y^2}\left[1(r\mathrm{cos}\varphi +x_0)^2(r\mathrm{sin}\varphi +y_0)^2\right]}`$ (15)
$`=`$ $`\mathrm{}n_0\stackrel{~}{R}_x^2\stackrel{~}{R}_y^2{\displaystyle 𝑑\varphi \frac{\dot{x}_0\mathrm{sin}\varphi +\dot{y}_0\mathrm{cos}\varphi }{\stackrel{~}{R}_x^2\mathrm{cos}^2\varphi +\stackrel{~}{R}_y^2\mathrm{sin}^2\varphi }\left(\frac{2}{3}B^3AB\right)},`$ (16)
where

$`A`$ $``$ $`1x_0^2y_0^2,`$ (17)
$`B`$ $``$ $`x_0\mathrm{cos}\varphi +y_0\mathrm{sin}\varphi ,`$ (18)
with $`x_0`$ and $`y_0`$ also dimensionless. As anticipated in the general discussion, $`𝒯`$ indeed depends on the first time derivatives of the parameters $`x_0`$ and $`y_0`$.

Lagrange’s equations now take the familiar form

$$𝒟_x\left(𝒯\right)=0,𝒟_y\left(𝒯\right)=0,$$
(19)
where

$$𝒟_x\frac{d}{dt}\frac{}{\dot{x}_0}\frac{}{x_0},$$
(20)
and similarly for $`𝒟_y`$. If this operator is applied to functions of the form appearing in Eq. (15), it is straightforward to verify that

$$𝒟_x\left[(\dot{x}_0\mathrm{sin}\varphi +\dot{y}_0\mathrm{cos}\varphi )f(x_0,y_0,\varphi )\right]=\dot{y}_0\left(\mathrm{cos}\varphi _{x_0}+\mathrm{sin}\varphi _{y_0}\right)f(x_0,y_0,\varphi ).$$
(21)
In our case,

$$\left(\mathrm{cos}\varphi _{x_0}+\mathrm{sin}\varphi _{y_0}\right)A=2B,\left(\mathrm{cos}\varphi _{x_0}+\mathrm{sin}\varphi _{y_0}\right)B=1.$$
(22)
Applying this operator to the time-dependent term (15) simplifies the integral considerably and yields

$`𝒟_x𝒯[\mathrm{\Psi }]`$ $`=`$ $`\dot{y}_0\mathrm{}n_0\stackrel{~}{R}_x^2\stackrel{~}{R}_y^2{\displaystyle 𝑑\varphi \frac{A}{\stackrel{~}{R}_x^2\mathrm{cos}^2\varphi +\stackrel{~}{R}_y^2\mathrm{sin}^2\varphi }}`$ (23)
$`=`$ $`\dot{y}_0\mathrm{\hspace{0.17em}2}\pi \mathrm{}n_0\stackrel{~}{R}_x\stackrel{~}{R}_y(1x_0^2y_0^2).`$ (24)
The first of Lagrange’s equations (19) then becomes

$$𝒟_x𝒯[\mathrm{\Psi }]=\frac{\mathrm{}[\mathrm{\Psi }]}{x_0},$$
(25)
where $`\mathrm{}`$ denotes the extra free energy associated with the presence of a vortex. This equation, along with a similar one for the $`y_0`$ dependence can be rewritten as

$$2\pi \mathrm{}n_0\stackrel{~}{R}_x\stackrel{~}{R}_y\widehat{z}\times \dot{𝐫}_0=\mathbf{}_0\mathrm{}(𝐫_0),$$
(26)
where $`𝐫_0=(x_0,y_0)`$ is the dimensionless coordinate vector of the displaced vortex. As is familiar in the context of two-dimensional vortices in superfluid helium , this dynamical equation can be interpreted as a balance between the Magnus force and the gradient of the free energy with respect to the coordinates of the vortex .

We can calculate the free-energy part with logarithmic accuracy  by substituting (9)-(11) into the energy functional (4):

$$\mathrm{}(𝐫_0,\mathrm{\Omega })=2\pi \stackrel{~}{\mu }\xi ^2n_0(1r_0^2)\left[\mathrm{ln}\left(\frac{\stackrel{~}{R}_{}}{\xi }\right)\frac{M\mathrm{\Omega }}{\mathrm{}}\frac{\stackrel{~}{R}_x^2\stackrel{~}{R}_y^2}{\stackrel{~}{R}_x^2+\stackrel{~}{R}_y^2}(1r_0^2)\right],$$
(27)
where $`\stackrel{~}{R}_{}^2=2\stackrel{~}{R}_x^2\stackrel{~}{R}_y^2/(\stackrel{~}{R}_x^2+\stackrel{~}{R}_y^2)`$ and $`\xi ^2=\mathrm{}^2/2M\stackrel{~}{\mu }`$. With these results, we can derive the equations of motion and obtain the angular frequency $`\omega _p(\mathrm{\Omega })`$ of the vortex precession. An easy calculation yields

$$\dot{x}_0=\omega _p(\mathrm{\Omega })y_0,\dot{y}_0=+\omega _p(\mathrm{\Omega })x_0,$$
(28)
where

$$\omega _p(\mathrm{\Omega })\omega _p^0\frac{2\stackrel{~}{R}_x\stackrel{~}{R}_y}{\stackrel{~}{R}_x^2+\stackrel{~}{R}_y^2}\mathrm{\Omega }$$
(29)
and

$$\omega _p^0\frac{\mathrm{}}{M\stackrel{~}{R}_x\stackrel{~}{R}_y}\mathrm{ln}\left(\frac{\stackrel{~}{R}_{}}{\xi }\right)\frac{1}{1r_0^2}$$
(30)
is the precession frequency in a nonrotating condensate of stretched dimensions $`\stackrel{~}{R}_x`$ and $`\stackrel{~}{R}_y`$. For a not-too-large angular velocity $`\mathrm{\Omega }\omega _{}`$, note that $`\omega _p^0\omega _p(0)`$. An earlier more intricate derivation of this result analyzed the dynamical motion of each element of the vortex core . If $`\omega _p(\mathrm{\Omega })`$ is positive, the vortex precesses in the positive (counterclockwise) sense, namely with the same sense as the circulating fluid around the vortex. Such behavior is indeed seen in the recent JILA experiments  for nonrotating condensates.

For an axisymmetric condensate with $`\stackrel{~}{R}_x=\stackrel{~}{R}_y=\stackrel{~}{R}_{}`$, the precession frequency $`\omega _p(0)`$ in the absence of rotation is shifted to $`\omega _p(\mathrm{\Omega })=\omega _p^0\mathrm{\Omega }`$ by the applied rotation (as expected from the transformation to the rotating frame ). For an asymmetric trap, however, the induced irrotational flow in Eq. (11) acts to screen the effect of the external rotation. In the extreme limit $`\stackrel{~}{R}_y\stackrel{~}{R}_x`$, for example, we have $`\omega _p(\mathrm{\Omega })\omega _p^0(2\stackrel{~}{R}_y/\stackrel{~}{R}_x)\mathrm{\Omega }`$, so that the precession then becomes nearly independent of $`\mathrm{\Omega }`$ for small $`\mathrm{\Omega }`$.

This behavior can be understood by examining the background fluid flow velocity in the rotating frame $`𝐯=(\mathrm{}/M)\mathbf{}^{}S_0\mathrm{\Omega }\widehat{z}\times 𝐫^{}`$ in the absence of the vortex ; for clarity, we now use primes to denote the original dimensional coordinates with, for example, $`x^{}=x\stackrel{~}{R}_x`$. In particular, $`v_x(x^{},y^{})=2\mathrm{\Omega }y^{}\stackrel{~}{R}_x^2/(\stackrel{~}{R}_x^2+\stackrel{~}{R}_y^2)`$. This induced velocity appears in the dimensional equation of motion as

$$\frac{\dot{x}_0^{}}{\stackrel{~}{R}_x}=\omega _p^0\frac{y_0^{}}{\stackrel{~}{R}_y}+\frac{v_x(x_0^{},y_0^{})}{\stackrel{~}{R}_x}=\omega _p^0\frac{y_0^{}}{\stackrel{~}{R}_y}+\frac{2\mathrm{\Omega }\stackrel{~}{R}_x}{\stackrel{~}{R}_x^2+\stackrel{~}{R}_y^2}y_0^{},$$
(31)
which readily reproduces the first of Eqs. (28), and similarly for the second equation.

## IV Three-dimensional condensate

It is not difficult to generalize these results to a three-dimensional TF condensate. For example, the explicitly time-dependent part of the Lagrangian now contains an additional integral over the axial coordinate $`z`$:

$`𝒯[\mathrm{\Psi }]`$ $`=`$ $`\mathrm{}n_0\stackrel{~}{R}_x^2\stackrel{~}{R}_y^2\stackrel{~}{R}_z{\displaystyle d^3r\frac{\dot{x}_0y+\dot{y}_0x}{\stackrel{~}{R}_x^2x^2+\stackrel{~}{R}_y^2y^2}\left[1(r\mathrm{cos}\varphi +x_0)^2(r\mathrm{sin}\varphi +y_0)^2z^2\right]}`$ (32)
$`=`$ $`2\mathrm{}n_0\stackrel{~}{R}_x^2\stackrel{~}{R}_y^2\stackrel{~}{R}_z{\displaystyle d^2r_0^{z_{\mathrm{max}}}𝑑z\frac{\dot{x}_0y+\dot{y}_0x}{\stackrel{~}{R}_x^2x^2+\stackrel{~}{R}_y^2y^2}\left[1(r\mathrm{cos}\varphi +x_0)^2(r\mathrm{sin}\varphi +y_0)^2z^2\right]}`$ (33)
$`=`$ $`{\displaystyle \frac{\mathrm{}n_0\stackrel{~}{R}_x^2\stackrel{~}{R}_y^2\stackrel{~}{R}_z}{6}}{\displaystyle 𝑑\varphi \frac{\dot{x}_0\mathrm{sin}\varphi +\dot{y}_0\mathrm{cos}\varphi }{\stackrel{~}{R}_x^2\mathrm{cos}^2\varphi +\stackrel{~}{R}_y^2\mathrm{sin}^2\varphi }}`$ (35)
$`\times \left\{5BA^{3/2}3B^3A^{1/2}+3\left(A+B^2\right)^2\left[{\displaystyle \frac{\pi }{2}}+\mathrm{arctan}\left({\displaystyle \frac{B}{A^{1/2}}}\right)\right]\right\},`$
where $`z_{\mathrm{max}}^2=1(r\mathrm{cos}\varphi +x_0)^2(r\mathrm{sin}\varphi +y_0)^2`$. Although this angular integral is difficult to evaluate directly, only the simpler quantity $`𝒟_x𝒯[\mathrm{\Psi }]`$ is needed, and a straightforward analysis yields

$`𝒟_x𝒯[\mathrm{\Psi }]=\dot{y}_0{\displaystyle \frac{8\pi \mathrm{}n_0}{3}}\stackrel{~}{R}_x\stackrel{~}{R}_y\stackrel{~}{R}_z(1r_0^2)^{3/2}.`$ (36)
In addition, the free energy is given by

$$\mathrm{}(𝐫_0,\mathrm{\Omega })=\frac{8\pi \stackrel{~}{\mu }\xi ^2n_0\stackrel{~}{R}_z}{3}(1r_0^2)^{3/2}\left[\mathrm{ln}\left(\frac{\stackrel{~}{R}_{}}{\xi }\right)\frac{4}{5}\frac{M\mathrm{\Omega }}{\mathrm{}}\frac{\stackrel{~}{R}_x^2\stackrel{~}{R}_y^2}{\stackrel{~}{R}_x^2+\stackrel{~}{R}_y^2}(1r_0^2)\right].$$
(37)
The dynamical equation again reduces to (26), and the only difference is that the angular frequency $`\omega _p^0`$ has a modified numerical coefficient

$$\omega _p(\mathrm{\Omega })\frac{3}{2}\frac{\mathrm{}}{M\stackrel{~}{R}_x\stackrel{~}{R}_y}\mathrm{ln}\left(\frac{\stackrel{~}{R}_{}}{\xi }\right)\frac{1}{1r_0^2}\frac{2\stackrel{~}{R}_x\stackrel{~}{R}_y}{\stackrel{~}{R}_x^2+\stackrel{~}{R}_y^2}\mathrm{\Omega };$$
(38)
this expression agrees with an earlier result based on the method of matched asymptotic expansions .

Note that our analysis for a three-dimensional condensate considers only a straight vortex line. In the present Thomas-Fermi limit of a large condensate, the axis of the vortex line must lie perpendicular to the condensate’s surface . Hence the approximation of a straight vortex is directly applicable to the central region of a disk-shaped condensate. For cigar-shaped condensates, in contrast, the bending of the vortex generally plays an essential role. Indeed, for sufficiently elongated condensates, a straight vortex on the central axis can become unstable with respect to bending deformations .

## V Discussion and conclusions

The Lagrangian functional (2) contains two parts, $`𝒯[\mathrm{\Psi }]`$ with explicit time dependence and $`[\mathrm{\Psi }]`$ involving the time-independent free energy in the rotating frame. If only the latter quantity is considered (which is the Hamiltonian in the rotating frame), the stability can be inferred directly by considering how $``$ changes for small lateral displacements of the vortex from the central position. Such methods have been used to analyze the onset of metastability in a rotating nonaxisymmetric TF condensate .

In contrast, the full Lagrangian allows a more complete dynamical description. In the present case, the time dependence of the vortex position $`𝐫_0`$ yields equations of motion (26), which readily give the precession frequency $`\omega _p(\mathrm{\Omega })`$ in a rotating condensate. As mentioned in Sec. II, the present description does not include an image vortex. Although such an image vortex is needed to satisfy the appropriate boundary condition for a condensate with uniform density in a rigid container , the present TF condensate density vanishes at the TF boundary so that the superfluid current automatically vanishes there. In addition, the image vortex has negligible effect on the energy when calculated with logarithmic accuracy. As emphasized in Refs. , the vortex is unstable (metastable) whenever the gradient of $`\mathrm{\Delta }(𝐫_0,\mathrm{\Omega })`$ with respect to $`𝐫_0`$ in Eq. (26) acts to move the vortex away from (back toward) the center of the trap.

The present Lagrangian approach gives a direct physical derivation of the precession frequency. Although the results describe only a single straight vortex displaced from the $`z`$ axis, the analysis is far simpler than the method of matched asymptotic expansions applied to the full Gross-Pitaevskii equation . For a one-component condensate, it fully justifies the intuitive dynamical equation (26) suggested earlier by McGee and Holland . It would be valuable to generalize the present approach to the more interesting and challenging case of two components, where additional restoring (buoyancy) forces have been proposed .

###### Acknowledgements.

This work has benefited from many helpful discussions with M. J. Holland and A. A. Svidzinsky. It has been supported in part by the NSF, Grant No. 99-71518.


text:
# Influence of Nd on the magnetic properties of 𝑁𝑑_{1-𝑥}𝐶𝑎_𝑥𝑀𝑛𝑂₃: an ESR study

## I Figure captions

Figure 1: temperature dependence of the DC magnetic susceptibility ($`\mathrm{}`$), ESR manganese susceptibility ($`\mathrm{\Delta }`$) normalized at room temperature and Nd susceptibility ($`\mathrm{}`$)calculated from ESR line shifts at $`0.33T`$ for $`Nd_{1x}Ca_xMnO_3`$, $`x=0.3`$ (Fig.1a) and $`x=0.5`$ (Fig.1b). The magnetic and charge order transitions have been determined by neutron diffraction at zero field.

Figure 2: temperature dependence of the Mn ($`\mathrm{}`$) and Nd ($`\mathrm{}`$) ordered ferromagnetic moment determined by neutron diffraction as well as the average ESR line shift ($``$) at $`95GHz`$ for $`Nd_{0.7}Ca_{0.3}MnO_3`$.

Figure 3: ESR absorption line position (corrected from the sample demagnetization field) as a function of temperature in $`Nd_{0.7}Ca_{0.3}MnO_3`$ taken at $`95GHz`$ (a) and in $`Nd_{0.5}Ca_{0.5}MnO_3`$ FM phase taken at $`475GHz`$ (b). the dotted line corresponds to $`g=1.99`$. Insert: ESR absorption spectrum at different temperatures.

Figure 4:  ESR absorption line position (corrected from the sample demagnetization field) as a function of temperature in $`Nd_{0.5}Ca_{0.5}MnO_3`$ taken at $`95GHz`$. The dotted line coresponds to $`g=1.99`$. Insert: ESR absorption spectrum at different temperatures.


text:
# Star Formation from Turbulent Fragmentation

## 1. Introduction

Stars are born in interstellar clouds of molecular hydrogen, and form by gravitational collapse of shock-compressed density fluctuations generated from the supersonic turbulence ubiquitously observed in molecular clouds (e.g. Elmegreen 1993, Padoan 1995, Klessen, Heitsch, & Mac Low 2000). Once a gas clump becomes gravitationally unstable, collapse sets in and the central density increases until a protostellar object forms which continues to grow in mass via accretion from the infalling envelope. As stars typically form in small aggregates or larger clusters (Lada 1992, Adams & Myers 2001) the interaction of protostellar cores and their competition for mass growth from their surroundings are important processes shaping the distribution of the final star properties.

Altogether star formation can be seen as a two-phase process: First, turbulent fragmentation leads to transient clumpy molecular cloud structure, with some of the density fluctuation exceeding the critical mass and density for gravitational contraction. This sets the stage for the second phase of star formation, the collapse of individual protostellar clumps building up central protostars. To fully understand the star formation process and all its phases in detail is one of the fundamental challenges in astronomy. Describing the first phase of star formation requires the three-dimensional dynamical modeling of molecular cloud formation, evolution and fragmentation. In this contribution, I report results from numerical simulations using SPH, a particle-based approach to solve the equations of hydrodynamics. The method and its application to (driven or decaying) interstellar turbulence is introduced in §2 and §3. I demonstrate how supersonic turbulence may influence the properties of forming (proto)stellar aggregates and clusters, and discuss the spatial distribution and timescale of star formation in §4 and the expected mass spectra for gas clumps and protostars in §5. Finally, in §6 a first attempt to combine all phases of star formation into a consistent numerical scheme is presented.

## 2. SPH – A Particle-Based Flexible Method to Hydrodynamics

To adequately describe turbulent fragmentation and the formation of protostellar cores, it is necessary to resolve the collapse of shock compressed regions over several orders of magnitude in density. Due to the stochastic nature of supersonic turbulence, it is not known in advance where and when local collapse occurs. Hence, SPH (smoothed particle hydrodynamics) is used to solve the equations of hydrodynamics. It is a Lagrangian method, where the fluid is represented by an ensemble of particles and thermodynamical observables are obtained by averaging over an appropriate subset of the SPH particles (Benz 1990). High density contrasts are resolved by simply increasing the particle concentration where needed. SPH can also be combined with the special-purpose hardware device GRAPE (Sugimoto et al. 1990, Ebisuzaki et al. 1993, Steinmetz 1996) allowing for calculations at supercomputer level on a normal workstation. The simulations presented here concentrate on subregions within a much larger cloud, therefore periodic boundary conditions are adopted (Klessen 1997). The high-density cores of collapsing gas clumps are substituted by ‘sink’ particles (Bate, Bonnell & Price 1995) while keeping track of mass and momentum. This allows for following the dynamical evolution of the system over many free-fall timescales.

## 3. Models of Driven Supersonic Turbulence

The large observed linewidths in molecular clouds imply the presence of supersonic velocity fields that carry enough energy to counterbalance gravity on global scales (Williams, Blitz, & McKee 2000). However, it is known that turbulent energy dissipates rapidly, roughly on the free-fall timescale (Mac Low et al. 1998, Stone, Ostriker, & Gammie 1998, Padoan & Nordlund 1999). Unlike previously thought, this is independent of the presence of magnetic fields. The fields are therefore taken as being dynamically unimportant in the current models, and are not included. To prevent or considerably postpone global collapse, turbulence must be continuously replenished. This is achieved by applying a Gaussian driving scheme, that inserts kinetic energy in a specified range of wavenumbers $`k`$. I select $`1k2`$, $`3k4`$, and $`8k8`$, corresponding to sources that act on large, intermediate, and small scales, respectively. The energy input at each timestep is adjusted to reach a constant level of the total kinetic energy sufficient to stabilize the cloud as a whole (Mac Low 1999). Furthermore one molecular cloud model is considered where turbulence is assumed to have decayed, leaving behind Gaussian density fluctuations which begin to contract on all scales. For details see Klessen (2001).

The models presented here are computed in normalized units. If scaled to mean densities of $`n(\mathrm{H}_2)=10^4`$cm<sup>-3</sup>, a value typical for star-forming molecular cloud regions (e.g. in $`\rho `$-Ophiuchus, see Motte, André, & Neri 1998) and a temperature of 11.4 K (i.e. a sound speed $`c_\mathrm{s}=0.2`$km/s), then the total mass contained in the computed volume is 413 M and the size of the cube is $`0.89`$pc. It contains 64 thermal Jeans masses.

## 4. Spatial Distribution and Timescale of Star Formation

Stars form from turbulent fragmentation of molecular cloud material. Supersonic turbulence that is strong enough to counterbalance gravity on global scales will usually provoke local collapse. Turbulence establishes a complex network of interacting shocks, where converging shockfronts generate clumps of high-density. This density enhancement can be large enough for the fluctuations to become gravitationally unstable and collapse, i.e. when the local Jeans length becomes smaller than the size of the fluctuation. However, the fluctuations in turbulent velocity fields are highly transient. The random flow that creates local density enhancements can disperse them again. For local collapse to actually result in the formation of stars, locally Jeans unstable, shock-generated, density fluctuations must collapse to sufficiently high densities on time scales shorter than the typical time interval between two successive shock passages. Only then are they able to ‘decouple’ from the ambient flow pattern and survive subsequent shock interactions. The shorter the time between shock passages, the less likely these fluctuations are to survive. Hence, the efficiency of protostellar core formation and the rate of continuing accretion onto collapsed cores depend strongly on the wavelength and strength of the driving source (see Klessen et al. 2000, Heitsch, Mac Low, & Klessen 2001).

The velocity field of long-wavelength turbulence is dominated by large-scale shocks which are very efficient in sweeping up molecular cloud material, thus creating massive coherent structures. When a coherent region reaches the critical density for gravitational collapse its mass typically exceeds the local Jeans limit by far. Inside the shock compressed region, the velocity dispersion is much smaller than in the ambient turbulent flow and the situation is similar to localized turbulent decay. Quickly a cluster of protostellar cores builds up. Both, decaying and large-scale turbulence lead to a clustered mode of star formation. The efficiency of turbulent fragmentation is reduced if the driving wavelength decreases. When energy is carried mainly on small spatial scales, the network of interacting shocks is very tightly knit, and protostellar cores form independently of each other at random locations throughout the cloud and at random times. Individual shock generated clumps have lower mass and the time interval between two shock passages through the same point in space is small. Hence, collapsing cores are easily destroyed again and star formation is inefficient. This scenario corresponds to the isolated mode of star formation. It needs to be pointed out that there is no fundamental dichotomy between the two modes of star formation, they rather define the extreme ends in the continuous spectrum of the properties of turbulent molecular cloud fragmentation.

This is visualized in Fig. 1. It compares the models of decayed, large-wavelength ($`k=1\mathrm{}2`$), intermediate ($`k=3\mathrm{}4`$), and small-scale turbulence ($`k=7\mathrm{}8`$). The density structure of the systems is depicted at $`t=0`$, which means for the model of decayed turbulence that the initial Gaussian density field is visible. For the driven models, $`t=0`$ corresponds to the phase of fully developed turbulence just before gravity is ‘switched on’. The lower panel describes the four models after the first protostellar cores have formed via turbulent fragmentation and have accreted 30% of the total mass. Time is measured in units of the global free-fall timescale. Dark dots indicate the location of dense collapsed core. As for decayed turbulence all spatial modes are unstable, the system quickly evolves into a network of intersecting filaments, where protostellar cores predominantly form. Similarly, also large-scale turbulence builds up a network of filaments, however, this time the large coherent structures are not caused by gravity, but instead are due to shock compression. Once gravity is included, it quickly dominates the evolution inside the dense regions and again a cluster of protostellar cores builds up. In the case of intermediate-wavelength turbulence, cores form in small aggregates, whereas small-scale turbulence leads to local collapse of individual objects randomly dispersed throughout the volume. Note the different times needed for 30% of the mass to be accumulated in dense cores. For small-scale turbulence star formation needs longest and speeds up with increasing wavelength, and as expected the rates for large-scale turbulence and locally decayed turbulence are comparable. This is also seen in Fig. 2, which shows the mass accretion history of individual protostellar cores for the three turbulent models, together with the distribution of core formation times. Star formation efficiency is high in the long-wavelength model, all cores form within two free-fall times, whereas in the short-wavelength model the efficiency is low and core formation continues for over 15 free-fall timescales (when the simulation was stopped).

## 5. Mass Spectra of Clumps and Protostellar Cores

The dominant parameter determining stellar evolution is the mass. It is therefore important to investigate the relation between the masses of molecular clumps, protostellar cores and the resulting stars. Figure 3 plots for the four models the mass distribution of all gas clumps, of the subset of Jeans critical clumps, and of collapsed cores. Four different evolutionary phases are shown, initially just when gravity is ‘switched on’, and after turbulent fragmentation has lead to the accumulation of $`M_{}5`$%, $`M_{}30`$% and $`M_{}60`$% of the total mass in dense cores.

In the initial, completely pre-stellar phase the clump mass spectrum is very steep (about Salpeter slope or less) at the high-mass end and gets shallower below $`M0.4M_\mathrm{J}`$ with slope $`1.5`$, when using a power-law fit. The spectrum strongly declines beyond the SPH resolution limit. Altogether, individual clumps are hardly more massive than a few $`M_\mathrm{J}`$.

Gravity strongly modifies the distribution of clump masses during the later evolution. As clumps merge and grow bigger, the mass spectrum becomes flatter and extends towards larger masses. Consequently the number of cores that exceed the Jeans limit grows, and local collapse sets in leading to the formation of dense cores. This is most evident in the Gaussian model of decayed turbulence, where the velocity field is entirely determined by gravitational contraction on all scales. The clump mass spectrum in intermediate phases of the evolution (i.e. when protostellar cores are forming but the overall gravitational potential is still dominated by non-accreted gas) exhibits a slope $`1.5`$ similar to the observed one. When the velocity field is dominated by strong (driven) turbulence, the effect of gravity on the clump mass spectrum is much weaker. It remains steep, close to or even below the Salpeter value. This is seen for small-wavelength turbulence. Here, the short interval between shock passages prohibits efficient merging and the build up of a large number of massive clumps. Only few clumps become Jeans unstable and collapse to form cores. These form independent of each other at random locations and times and typically do not interact. Increasing the driving wavelength leads to more coherent and rapid core formation, which also results in a larger number of cores. It is maximum in the case of pure gravitational contraction (decayed turbulence).

Long-wavelength turbulence or turbulent decay leads to a core mass spectrum that is well approximated by a log-normal. It roughly peaks at the average thermal Jeans mass $`M_\mathrm{J}`$ of the system and is comparable in width with the observed IMF (see Klessen & Burkert 2000, 2001). The log-normal shape of the mass distribution may be explained by invoking the central limit theorem (e.g. Zinnecker 1984), as protostellar cores form and evolve through a sequence of highly stochastic events (resulting from supersonic turbulence and/or competitive accretion). To find the mass peak at $`M_\mathrm{J}`$ may be somewhat surprising given the fact that the local Jeans mass strongly varies between different clumps. In a statistical sense the system retains knowledge of its mean properties. The total width of the core distribution is about two orders of magnitude in mass and is approximately the same for all four models. However, the spectrum for intermediate and short-wavelength turbulence, i.e. for isolated core formation, is too flat (or equivalently too wide) to be comparable to the observed IMF. This is in agreement with the hypothesis that most stars form in aggregates or clusters.

## 6. Towards an Integral Description of Star Formation

Three-dimensional calculations of molecular cloud dynamics can only describe the first phases of the star formation process: turbulent molecular cloud fragmentation and formation of collapsing protostellar cores. These calculation fail to resolve the build-up of the protostar itself. Understanding the temporal evolution of its structural parameters and the determination of observables like luminosity and effective temperature as the star approaches the main sequence requires detailed protostellar evolutionary calculations. These calculations are important, because they give ages, masses and radii of young stars when brightness, distance and effective temperature are known. The so determined ages constitute the only practical ‘clock’ for tracing the history of star-formation regions and for studying the evolution of circumstellar disks and planet formation.

However, almost all classical PMS calculations start after the main protostellar accretion phase ended. Formation and mass growth of the protostar are usually neglected, and its internal thermal structure is taken to be fully convective (e.g. D’Antona & Mazzitelli 1994). However, it has been shown recently that protostellar collapse does not lead to a fully convective stellar structure when calculated consistently, taking into account radiation-hydrodynamics, time-dependent convection and nuclear burning processes (Wuchterl & Tscharnuter 2001).

A fully consistent treatment of formation and evolution of protostars requires coupling the 3D cloud fragmentation models with the 1D dynamical description of protostellar collapse. This can be achieved by taking the physical parameters and the time dependent mass accretion rates onto individual protostellar cores from the 3D model as boundary values for the detailed 1D calculation (see Wuchterl & Klessen 2001). For the proto-Sun, i.e. a cloud fragment with $`1`$M, the evolutionary path in the Hertzsprung-Russell (HR) diagram is illustrated in Fig. 4.

These calculations demonstrate that luminosity and temperature of protostars not only depend on mass and age, but also on the accretion rate in the main accretion phase. Therefore, the position of a protostar on the HR diagram is sensitive to the cluster environment, because stochastic core interaction and competition strongly influence the mass accretion rate: PMS tracks are not unique, and mass (or similarly age) of individual protostars from using the HR diagram can only be determined in a statistical sense.

Furthermore, the new dynamical models predict considerable differences to the currently established classical PMS tracks: Instead of being fully convective, the structure of a $`1`$M protostar is homologous to the present day Sun. Already an age of 1 million years, it has a radiative core and a convective envelope. For this reason, the dynamical PMS calculacions predict that the proto-Sun was $`500`$K hotter and twice as luminous as in the classical hydrostatic model. Altogether, conclusive observational tests are necessary to determine the validity of the numerical models. It is essential in this context to obtain independent mass and/or age determinations, e.g. from PMS binary stars or detailed modeling of protostellar disks (see Simon, Dutrey, & Guilloteau 2000, or Woitas, Köhler, & Leinert 2001 for first estimates).

### Acknowledgments.

I thank Günther Wuchterl, Mordecai-Mark Mac Low, Fabian Heitsch, and Peter Bodenheimer for many stimulating discussions and fruitful collaboration.

## References

Adams, F. C., Myers, P. C., 2001, ApJ, 553, 744

Benz, W., 1990, in The Numerical Modeling of Nonlinear Stellar Pulsations, ed. J. R. Buchler (Dordrecht: Kluwer), 269

Bate, M. R., Bonnell, I. A., Price, N. M., 1995, MNRAS, 277, 362

D’Antona, F., Mazzitelli, I., 1994, ApJS, 90, 467

Ebisuzaki, T., Makino, J., Fukushige, T., Taiji, M., Sugimoto, D., Ito, T., Okumura, S. K., 1993, PASJ, 45, 269

Elmegreen, B. G., 1993, ApJ, 419, L29

Heitsch, F., Mac Low, M.-M., Klessen, R. S., 2001, ApJ, 547, 280

Klessen, R. S., 1997, MNRAS, 292, 11

Klessen, R. S., 2001, ApJ, 556, 837

Klessen, R. S., Burkert, A., 2000, ApJS, 128, 287

Klessen, R. S., Burkert, A., 2001, ApJ, 549, 386

Klessen, R. S., Heitsch, F., Mac Low, M.-M., 2000, ApJ, 535, 887

Lada, E., 1992, ApJ, 393, L25

Mac Low, M.-M., 1999, ApJ, 524, 169

Mac Low, M.-M., Klessen, R. S., Burkert, A., Smith, M. D., 1998, Phys.Rev.Lett, 80, 2754

Mizuno, A., Onishi, T., Yonekura, Y., Nagahama, T., Ogawa, H., Fukui, Y., 1995, ApJ, 445, L161

Motte, F., André, P., Neri, R., 1998, A&A, 336, 150

Padoan, P., 1995, MNRAS, 277, 337

Padoan, P., Nordlund, Å., 1999, ApJ, 526, 279

Steinmetz, M., 1996, MNRAS, 278, 1005

Simon, M., Dutrey, A., Guilloteau, S., 2000, ApJ, 545, 1034

Stone, J. M., Ostriker, E. C., Gammie, C. F., 1998, ApJ, 508, L99

Sugimoto, D., Chikada, Y., Makino, J., Ito, T., Ebisuzaki, T., Umemura, M., 1990, Nature, 345, 33

Williams, J. P., Blitz, L., McKee, C. F., 2000, in Protostars and Planets IV, eds. V. Mannings, A. Boss, & S. Russell (Tucson: Univ. of Arizona Press), 97

Woitas, J., Köhler, R., Leinert, C., 2001, A&A, 369, 249

Wuchterl, G., Tscharnuter, W. M., 2001, A&A, submitted

Wuchterl, G., Klessen, R. S., 2001, ApJ, 560, in press (astro-ph/0109051)

Zinnecker, H., 1984, MNRAS, 210, 43


text:
# Detection of spiral structure of the quiescent accretion disk of IP Pegasi

## References

Steeghs D., Harlaftis E.T., Horne K., 1997, MNRAS, 290, 28

Wolf S., Barwig H., Bobinger A. et al., 1998, A&A, 332, 984