text:
# Magnetism in Atomic-Sized Palladium Contacts and Nanowires

## Abstract

We have investigated Pd nanowires theoretically, and found that, unlike either metallic or free atomic Pd, they exhibit Hund’s rule magnetism. In long, monoatomic wires, we find a spin moment of 0.7 $`\mu _B`$ per atom, whereas for short, monoatomic wires between bulk leads, the predicted moment is about 0.3 $`\mu _B`$ per wire atom. In contrast, a coaxial (6,1) wire was found to be nonmagnetic. The origin of the wire magnetism is analyzed.

Magnetism at the nanoscale is an exciting emerging research field, of both basic and applied relevance. In future technology, when the electronic components will have become so tiny that quantum size effects will determine their functionality, an understanding of nanomagnetic phenomena will be crucial; yet, relatively little is currently understood about how magnetism arises and how it affects the properties of metals at the nanoscale.

Systems of special interest are nanowires and atomic-sized nanocontacts. The low dimensionality of such systems causes specific physical phenomena to appear, for example quantized ballistic conductance wees1988  and helical wire geometries gulseren1998 . These phenomena interplay with the possible presence of magnetism in the nanosystem, especially of a genuine Hund’s rule magnetic order parameter. Here, we report theoretical studies of emerging magnetism in Pd nanowires in various geometries, see Fig. 1. Monowires of Pd, i.e., wires consisting of a single line of atoms, have recently been observed by Rodrigues et al.rodrigues . We find magnetic moments as high as 0.7 $`\mu _\mathrm{B}`$ per atom for infinitely long, straight, monatomic wires. Even short, three-atom Pd wire chains between bulk-like leads, are predicted to possess a magnetic moment around 0.3 $`\mu _\mathrm{B}`$ per wire atom, whereas the bulk leads remain nonmagnetic. In contrast, thicker coaxial (6,1) wires were found to be nonmagnetic.

Of course, thermal fluctuations, very large in a nanosystem, will generally act to destroy static magnetic order in the absence of an external field. There are nonetheless two different fluctuation regimes: slow and fast. Slow fluctuations transform a nanomagnet to a superparamagnetic state, where magnetization fluctuates on a long time scale between equivalent magnetic valleys, separated, e.g., by anisotropy-induced energy barriers. If the barriers are sufficiently large, the nanosystem spends most of the time in a single magnetic valley, and will for many practical purposes behave as magnetic. We may under these circumstances be allowed to neglect fluctuations altogether, and to approximate (as we will do here) the calculated properties of the superparamagnetic nanosystem with those of a statically magnetized one. Experimentally, evidence of one-dimensonal (1D) superparamagnetism with fluctuations sufficiently slow on the time scale of the probe has been reported in Co atomic chains deposited at Cu surface steps gambardella2002 . At the opposite extreme — a situation reached for example at high temperatures — the energy barriers are so readily overcome that the magnetic state is totally washed away by fast fluctuations, leading to a conventional paramagnetic state, in which all magnetic exchange splittings of the electronic bands vanish through motional narrowing. A complete description of this high entropy state is beyond scope, but a conventional $`T=0`$ nonmagnetic, singlet solution of the Kohn-Sham electronic structure equations can be used in its place, at least as a crude approximation.

The present density functional calculations dft  were performed with an all-electron full-potential linear muffin-tin orbital (FP-LMTO) basis set wills  together with a generalized gradient approximation (GGA) gga  to the exchange-correlation functional. A simpler local density functional lda  was also tested, giving results very similar the GGA. As a double check, some of the calculations were repeated using the linear augmented plane-wave code WIEN97 wien97 . None of these calculations assume any shape approximation of the potential or wave functions. We performed both scalar relativistic calculations, and calculations including the spin-orbit coupling as well as the scalar-relativistic terms. In the calculations involving the spin-orbit interaction, the spin axis was chosen to be aligned along the wire direction.

The calculations were performed with inherently three-dimensional codes. Thus, the infinitely long, straight, isolated monowires, as well as the (6,1) coaxial wires, were simulated by regular arrays of well spaced nanowires. The short wire consisted of three atoms in a straight line attaching to thick planar slabs of close-packed bulk Pd. The bond length in the short wire was chosen to be 2.7 Å, as suggested by transmission electron microscope images of Pd monowires rodrigues , and corresponds to a somewhat stretched wire. Convergence of the magnetic moment was checked with respect to k-point mesh density, Fourier mesh density, tail energies, wire-wire vacuum distance and bulk thickness. In the calculations, spin-orbit coupling was seen to have only a very minor effect on the results, in contrast to the situation for $`5d`$ metalsdelin\_wires .

Fig. 2 shows the magnetic moment per atom for the (infinitely) long Pd monowire, as a function of the bond length along the wire. The total magnetic moment rapidly reaches 0.7 $`\mu _B`$ at a bond length around 2.3 Å, and retains this magnetic moment over a long region. The steep rise at the onset of magnetism is almost entirely due to $`4d`$ polarization, which reaches a maximum of $``$0.5 $`\mu _B`$ already at 2.6 Å and then decreases monotonically. The remainder of the magnetic moment has $`5s`$ character. At around 3.5 Å, the magnetic moment disappears completely, where the wire undergoes a metal-insulator transition, with the opening of a $`ds`$ gap, foreshadowing that of the noninteracting atoms.

For the long monatomic wire, the energy gain per atom due to spin polarization is around 25 meV at the equilibrium bond length 2.56 Å. For comparison, we mention that in bulk Ni the corresponding energy gain is around 40 meV, indicating that the Pd wire magnetism might actually exist not only at ultra-low temperatures. Antiferromagnetic Pd monowire configurations were also tested, and found to be energetically unstable compared to the ferromagnetic configuration. Our result differs from that of Bahn et al.bahn2001  who found no magnetism in pseudopotential calculations for Pd monowires. It is possible that the disagreement could arise in this very borderline case due to the different methods used, in which case we would tend to trust our all-electron approach better.

Why does the Pd monowire magnetize? In the bulk, and also at surfaces, the Pd $`4d`$ band is too wide to provoke spin-polarization. In the atom, on the other hand, Pd forms a singlet with a completely filled $`4d`$ shell. Thus, the limiting cases are all nonmagnetic and the strong Hund’s rule magnetism in the wire appears rather unexpected. However, the borderline case of Pd magnetism is demonstrated by the fact that certain clusters are predicted to have a magnetic ground state landman2001 . From the atomic perspective, as the distance between Pd atoms decreases, the increased interatomic hybridization causes the $`4d`$ band to become partially unfilled due to $`ds`$ transfer, and thus symmetry breaking through spin-polarization becomes possible. Hund’s rules make it reasonable to assume that spin polarization will also be energetically favorable. From the bulk perspective, the reduction of the number of nearest neighbors in the wire compared to the bulk causes a narrowing of the $`4d`$ band, and the band width may become sufficiently small that the gain in exchange energy due to spin polarization is larger than the increase in kinetic energy. Thus, using this $`4d`$ band-narrowing argument, we can rationalize the existence of a magnetic state in the wire.

However, we have reason to believe that this explanation of the large magnetic moment in long Pd monowires is too simple, and that $`4d`$ band-narrowing is not the only mechanism at work, but that the one-dimensionality of the system is crucial. We will illustrate this point by analyzing the band structure of the long Pd monowires in some detail.

Fig. 3 shows the band structure of the Pd monowire for several different bond lengths, as it goes from insulator, to ferromagnet, to paramagnet, for decreasing bond length. There are six orbitals to consider, two $`s+d_{z^2}`$ orbitals, two $`(d_{xz},d_{yz})`$ orbitals, and two $`(d_{xy},d_{x^2y^2})`$ orbitals, as illustrated in Fig. 4. The two $`s+d_{z^2}`$ orbitals give rise to the two bands with highest dispersion. Of these, the band with mostly $`s`$ character has an electron band edge close to the Fermi level at $`\mathrm{\Gamma }`$, and the $`d`$-character dominated one has a hole band edge close to the Fermi level at the zone boundary A. Next, we have the two $`(d_{xz},d_{yz})`$ orbitals, which are degenerate and have a band edge close to the Fermi level at $`\mathrm{\Gamma }`$. Finally, the two $`(d_{xy},d_{x^2y^2})`$ orbitals form a flat band situated between 0.5 and 1 eV below the Fermi level.

The $`(d_{xy},d_{x^2y^2})`$ band lies far below the Fermi level with both spins completely occupied, and gives no contribution to magnetism. The spin moment is related to the $`s+d_{z^2}`$ and $`(d_{xz},d_{yz})`$ bands, which have a high dispersion, and display one-dimensional band edges close to the Fermi level at $`\mathrm{\Gamma }`$ and A. In the magnetic regime, i.e., from 2.3 Å to 3.5 Å, these three band edges all are nearly degenerate in energy and close to the Fermi level. This accidental feature of the long Pd monowire band structure dramatically increases the density of states, with divergent van Hove singularities close to the Fermi level. A spin splitting of these band edges can thus lower the total electron band energy, somewhat analogous to a band Jahn Teller, or a Peierls instability effect.

The fact that the Stoner stability criterion against magnetism is necessarily violated leading to a magnetized nanowire when a band edge approaches to the Fermi level is universal, as pointed out earlier by Zabala et al. zabala1998 . However, whereas for a simple $`sp`$ metal the ensuing 1D magnetism will imply very small moments which will readily be washed away by fluctuations, here the situation is different. The moment is large, and we expect a robust superparamagnetic state for temperatures not too high, transformable to a genuine magnetic state under an external field. The 1D singularities — in scattering terms the relevant band edge singularities correspond to standing waves, here of $`s+d_{z^2}`$ and $`(d_{xz},d_{yz})`$ electrons and holes — provide here merely the mechanism that triggers awake a strong Hund’s rule magnetic moment.

If the Pd monowire magnetism is to the most part due to a 1D scattering phenomenon, it should disappear or at least become much smaller when the number of nearest neighbors becomes large. To test this assertion, we, as already mentioned, also performed calculations for a two-layer coaxial (6,1) wire. As expected, it was found to be nonmagnetic, supporting the view that the magnetism in the Pd monowire is, to a large extent, a 1D phenomenon. This 1D magnetism in Pd proves nevertheless to be remarkably robust, as it survives even down to three-atom chains suspended between nonmagnetic leads.

The magnetic moment in the short wire suspended between bulk Pd leads is substantially smaller than the moment of the long monowire moment; 0.3 $`\mu _B`$ per wire atom compared to 0.7 $`\mu _B`$. Due to the high complexity of the band structure for the short wire (not shown), it is not possible to single out band edges as responsible for triggering the Hund’s rule magnetism in the short-wire case. However, we find that the bulk leads affect the electronic structure in the short wire resulting in a wider $`4d`$ band than in the long wire, which in turn reduces the magnitude of the magnetic moment.

In summary, we predict that ultimately thin Pd nanowires, both short and long, exhibit a spin-polarized ground state. The resulting superparamagnetic state of the nanowire should show up in the ballistic conductance in the form of a strong and unusual magnetic field and temperature dependence as well as spin-polarization of the current through the wire. Rodrigues et al.rodrigues  recently measured the charge conductance of Pd nanocontacts and found features in the conductance histogram above as well as below one conductance quantum $`G_0=2e^2/h`$. The existence of features below $`G_0`$ is intriguing since it seems to suggest that the $`s`$-dominated band crossing the Fermi level is spin-split around the Fermi level. A speculative explanation could entail some kind of spin reversal amid the nanowire smogunov2002 . More theory work will be needed to address their data, explicitly including such elements as tips, spin structures, and temperature as well as their effects on the system’s conductance. We are currently working to address this issue.

###### Acknowledgements.

A.D. acknowledges financial support from the European Commission through contract no. HPMF-CT-2000-00827 Marie Curie fellowship, STINT (The Swedish Foundation for International Cooperation in Research and Higher Education), and NFR (Naturvetenskapliga forskingsrådet). Work at SISSA was also sponsored through TMR FULPROP, MUIR (COFIN and FIRB RBAU01LX5H) and by INFM/F. J. M. Wills is acknowledged for letting us use his FP-LMTO code. We are also grateful to D. Ugarte for sharing with us the unpublished results of Ref. rodrigues, .


format: text

text:
# KH 15D: a star eclipsed by a large-scale dusty vortex?

## 1  Introduction

KH 15D is a pre-main sequence (PMS) K7 star located in the young star cluster, NGC 2264. Its large photometric variations are periodic and possibly interpreted in terms of eclipses by an opaque feature orbiting inside a circumstellar disk seen nearly edge-on (Herbst et al., 2002). It is likely connected to the UX Ori star family, characterized by large photometric and polarimetric variabilities (eg. Natta et al., 1999), but sets apart from them with two important differences: no infrared excess in the disk and no significant color changes during the eclipses, both indicating particle sizes much larger than optical wavelengths. As the $`VRI`$ colours of the star are consistent with its spectral type (Hamilton et al., 2001), the selective absorption is $`A_\mathrm{V}<0.2`$, a value indeed in agreement with the foreground interstellar reddening, $`E_{\mathrm{B}\mathrm{V}}0.060.07`$, of stars in NGC 2264 (Pérez et al., 1987; Sung et al., 1997). On the other hand the grey absorption, $`A_0`$, during the out-of-eclipse observations can be estimated indirectly using a well known relation (Barnes et al., 1978) between the visual brightness parameter and the intrinsic colour $`(VR)_0`$. Using the absolute visual magnitude deduced from the distance modulus of the cluster, this relation provides a value of the stellar radius consistent with that derived by Hamilton et al. (2001) and with the PMS status only if the grey absorption is $`0.60.7`$ magnitude.

The five eclipses observed since the discovery of KH 15D (Kearns and Herbst, 1998) provide us with the main characteristics of the occultation: a relative depth $`0.95`$, an approximate duration of 18 days and a period of 48.36 days. The absorbing cloud is deduced to orbit at $`0.2`$ AU from the star and to extend over 1/3 of an orbit (Hamilton et al., 2001). Such an occultation completely differs from a standard planetary transit and also displays other important peculiarities: (i) ingress and egress are remarkably steep, indicating an occulting clump with a very sharp edge; (ii) the minimum is not a flat land but contains a central bump with small scale features and a reverse peak close to mid-eclipse; (iii) when comparing two successive occultations, the reverse peak seems to shift from one side of the mid-eclipse to the other; (iv) since the first series of observations, the occultation has clearly widened and deepened in secular fashion.

The various models suggested by Herbst et al. (2002) to explain these observations involve a companion, planet or brown dwarf, able to $`\mathrm{`}\mathrm{`}`$shepherd” gravitationally an extended swarm of solid particles (warped density wave or ring arc) which is responsible for the eclipses. Possibly, these models could account for a number of the observed characteristics but, so far, none of them has been developped far enough to reproduce the observed light-curve and no observation has confirmed the presence of a $`\mathrm{`}\mathrm{`}`$shepherd”.

In this letter, we propose a new model in which KH 15D is periodically occulted by a dusty anticyclonic vortex persisting for a long time in a gas disk surrounding the star. The particles are confined inside the vortex by the gas friction and no gravitational assistance is needed to keep them together into a clump. The swarm of the trapped particles is organized by the gas dynamics with an absorption profile that pretty well explains most properties of the observed light-curve. The gaseous vortex associated to this dust swarm covers a wide azimuthal extent (see Fig.1) which seems quite puzzling to explain but looks like those produced in recent numerical simulations (Klahr and Bodenheimer, 2003). Our model is based on a scenario proposed by Barge and Sommeria (1994,1995; thereafter BS94,BS95) in which the formation of planets begins inside persistent gaseous vortices; it relies on a number of recent works on the formation and evolution of vortices in circumstellar disks (Li et al., 2001; Klahr and Bodenheimer, 2003).

## 2  The proposed model

Our model starts from a circumstellar disk of gas and solid particles in which the mechanisms of planetary formation are still at work. The gas is assumed to spread into a flared disk following a simple hydrostatic equilibrium in the vertical direction with a scale height $`H=C_\mathrm{S}/\mathrm{\Omega }`$, where $`C_\mathrm{S}`$ is the sound speed and $`\mathrm{\Omega }`$ is the Keplerian orbital frequency. Protoplanetary disks are likely turbulent during a span of their live and could host large scale and long-lived vortices growing from turbulence by analogy with what happens in two-dimensional fluid dynamics, in which organized structures are known to emerge from random turbulence in rotating shear flows. Vortices spinning like the shear flow are robust and merge one another while those with the opposite sign are laminated by the shear (BS95 and references therein). However, such vortices may also result from some specific instabilities as for example: (i) the Rossby Wave instability investigated by Li et al. (2001) which requires a strong local maximum in entropy and density; (ii) a global baroclinic instability described recently by Klahr and Bodenheimer (2003) and which arises from a natural radial stratification of the gas flow.

A number of numerical simulations show that large scale vortices can grow and survive for many rotation periods and look like either elongated structures stretched by compressibility effects (Godon and Livio, 2000) or a single vortex dominating the whole disk (Klahr and Bodenheimer, 2003). Highly elongated vortices are found also in the inner regions of MHD accretion disks (Tagger and Pellat 1998, and private communication) or as exact solutions of the incompressible Euler equation with different aspect ratios (Chavanis, 2000). On the other hand, the radial extent of vortices cannot exceed the thickness of the disk ($`R<H`$) since the velocity of the vortex, $`R\mathrm{\Omega }`$, must be less than the sound speed $`C_\mathrm{S}=H\mathrm{\Omega }`$ (BS95). Larger vortices would be destroyed by energy losses due to sound waves and density waves. The solid particles embedded in the gas disk are submitted to a friction drag depending on the mean-free-path $`\lambda `$ of the gas molecules relative to the particle size (or radius) $`s`$. In the inner disk region ($`r0.2`$AU) and for particles larger than a millimeter, the mean-free-path is less than size and the drag, caused by particle wake, reaches the Stokes regime (in the limit of particle Reynolds numbers $`Re_\mathrm{p}<1`$) with a stopping-time (Chavanis, 2000):

$$T_s=\frac{8s^2\rho _0}{9\sigma _{\mathrm{gas}}\mathrm{\Omega }\lambda },$$
(1)
where $`\sigma _{\mathrm{gas}}`$ is the gas surface density and $`\rho _02\mathrm{g}/\mathrm{cm}^3`$ is the density of material particles are made of. In contrast, the Epstein regime is reached for smaller particle size $`s<9\lambda /4`$. The particles are only submitted to star attraction and gas drag; their dynamical evolution depends on the single non-dimensioned friction parameter $`\tau _\mathrm{s}=\mathrm{\Omega }T_\mathrm{s}`$: (i) the lightest particles ($`\tau _\mathrm{s}1`$) come at rest rapidly with the gas and are driven by the flow; (ii) the heaviest particles ($`\tau _\mathrm{s}1`$) are nearly unaffected by the gas motion and keep a quasi-Keplerian motion.

The capture of particles by a gaseous vortex has been explored first by BS94 and BS95 using a simple model in which the velocity field is made up of concentric epicycles inside the vortex, while it matches a Keplerian flow at large distances. They found that such an anticyclonic vortex can capture and concentrate dust particles very efficiently: (i) light particles penetrate into the vortex and stop on streamlines close to the edge where they slowly shift toward the core; (iii) optimal particles, with $`\tau _\mathrm{s}1`$, sink deeply into the vortex and reach core streamlines. These results were confirmed for other velocity fields (Chavanis, 2000; De la Fuente et al., 2002). The vortex induces a segregation of the trapped particles (or a size sorting if the particles have the same composition) following the value of their friction parameter; this is illustrated in a number of numerical simulations (De la Fuente et al., 2002).

The ultimate reason for this dynamical behaviour lies in the sign of the vortex rotation. Indeed, in a reference frame rotating with the vortex center, the Coriolis force can overcome the centrifugal force and pushes the particles toward the core if the vortex is an anticyclone, whereas both forces are conspiring to eject the particles for a cyclone. This capture-in-vortex mechanism is a very efficient one and results in strong density enhancements inside the vortex, by at least two order of magnitude in $`200`$ rotation periods. The capture rate is estimated under the assumption that the particles are continuously renewed near the vortex orbit due to the inward drift under the systematic headwind drag (Weidenschilling 1977, BS95). One obvious and important consequence of this density enhancement is that, inside a vortex, particle growth is made easier and will couple with confinement and segregation. On the other hand, dust is depleted from the region inside the vortex orbit as the particles either are feeding the vortex or are falling to the star under the systematic drift.

Inside the vortex, the trapped particles are also submitted to a background small scale turbulence which makes them diffuse and tends to reduce their global concentration. Chavanis (2000) investigated this question in terms of a diffusion equation in an idealized circular vortex and derived a time dependant solution for the surface density inside the vortex:

$$\sigma _\mathrm{d}\frac{1}{l_\mathrm{d}^2(1k^2)}\mathrm{exp}\left[\frac{(rkr_0)^2}{l_\mathrm{d}^2(1k^2)}\right],$$
(2)
the initial state being a delta function centered at $`r_0`$; $`l_\mathrm{d}`$ is the diffusion scale length and $`k=\mathrm{exp}(t/T_{\mathrm{capt}})`$ in which $`T_{\mathrm{capt}}`$ is the characteristic time for a particle to reach the center of the vortex. In the case of light particles and very elongated vortices $`T_{\mathrm{capt}}2q/(3\mathrm{\Omega }\tau _\mathrm{s})`$ and $`l_d\sqrt{\alpha _\mathrm{v}/q\tau _\mathrm{s}}R_\mathrm{v}`$ where $`R_\mathrm{v}`$ is the vortex radius and $`q`$ its aspect ratio; $`\alpha _\mathrm{v}`$ is the non-dimensioned parameter measuring the small scale turbulence efficiency inside the vortex.

In order to estimate the size of the trapped particles, we will choose a standard model of nebula, the minimum mass solar nebula, in which the surface densities (both for gas and particles) and the temperature are the decreasing power laws $`r^{3/2}`$ and $`r^{1/2}`$, respectively; at 1 AU the densities are set to 1700 $`\mathrm{g}\mathrm{cm}^2`$ for the gas and 20 $`\mathrm{g}\mathrm{cm}^2`$ for the particles, whereas the temperature is assumed to be 280 K. At 0.2 AU from the star and with our numerical values, the optimal size for particle capture is $`s_{\mathrm{opt}}31\sqrt{\tau _\mathrm{s}}(r/1\mathrm{A}\mathrm{U})^{5/8}11`$ cm, i.e. pebble size. Lighter non optimal particles with $`\tau _\mathrm{s}<1`$ can be also captured by the vortex but remain trapped at the periphery of the vortex. These particles remain in the Stokes regime, at 0.2 AU from the star, as far as their size is larger than the critical size $`s_\mathrm{c}=9\lambda /40.3`$ mm.

## 3  Fitting the light-curve

We assume the disk is seen nearly edge on, under an inclination $`i`$ less than the flaring angle of the gas disk. $`i`$ must be small enough for the line of sight to cross the vortex in its vertical extent, but also large enough to avoid the prohibitive optical depths near the mid-plane layers. Of course, this inclination has to be consistent with the out-of-eclipse grey extinction $`A_00.7`$ magnitude deduced in Sect.1. In its motion around the star, the vortex periodically crosses the line of sight and the optical depth $`\tau `$ is a function of time varying from $`A_0`$ (out of eclipse) to $`A_0+3.2`$ during the eclipses. Assuming the absorption is due to spherical particles with sizes larger than a millimeter, the opacity $`\kappa `$ reduces to $`\pi s^2n`$ and the optical depth $`\tau =\kappa 𝑑l`$ is the familiar dust column density, where $`n`$ is the local number density along the line of sight.

During an occultation $`\tau `$ is found to increase inside the vortex, from center to edge, then to fall down very steeply to reach the out-of-vortex level. This $`\mathrm{`}\mathrm{`}`$opacity-curve” has a peculiar shape that our model can easily reproduce.

First, one can guess that a circular blob of matter in which the opacity $`\kappa (r)`$ is radially increasing can conveniently mimic the optical depth profile. Then, we noticed that a simple linear dependance, $`\kappa 1+ar`$, with a cut off at the outer boundary permits to fit pretty well the data (Fig. 2). In our fit the optical depth has been computed following two steps: (i) an integration of $`\kappa `$ along a path crossing an equivalent circular vortex (ECV) of radius $`H`$ and located at a distance $`r_{\mathrm{orb}}`$ from the star; (ii) a circular anamorphosis by a factor $`q17`$ (the result of an integration along a path at angle $`u`$ is reported to the direction of observation at angle $`q\times u`$). Figure 1 illustrates the ECV model with a sketch of the integration path and the density contours of the vortex based on the linear opacity law and the appropriate anamorphosis.

This absorption profile corresponds to a shell-like distribution of the solid particles inside the vortex, reminiscent of the spherical structures observed in Hii regions. In our model such a density profile originates in the dynamics and the segregation of particles trapped in a vortex as described above.

Surprisingly, a simple cut off at the vortex outer boundary can easily mimic the very steep fall during ingress or egress, consistent with an occultation by a sharp edge as proposed by Herbst et al. (2002). In fact, a cut off cannot receive true physical justification and the resulting profile cannot connect smoothly the out-of-eclipse level. A more realistic fit of ingress and egress directly arises from our model. Indeed, light sub-optimal particles are preferentially captured in the outer part of the vortex where they accumulate and slowly diffuse toward the core. The fit is realized thanks to the Gaussian density profile presented in Eq. 2 with the following values of the parameters: the center of the Gaussian is close to the vortex boundary (at a distance from the core $`R_\mathrm{v}=qH`$); the $`1/e`$ half-width of the Gaussian is given by $`l_\mathrm{d}\sqrt{(1k^2)}0.07R_\mathrm{v}`$. This is possible using the following assumptions: (1) the particles are renewed by capture every $`10`$ rotations of the vortex and have a friction parameter $`\tau _\mathrm{s}0.01`$, consistent with a capture at the vortex periphery and a slow drift toward the core in a time scale $`T_{\mathrm{capt}}2q/(3\mathrm{\Omega }\tau _\mathrm{s})`$; (2) the parameter measuring the turbulence inside the vortex is $`\alpha _\mathrm{v}0.01`$. As a result, the correponding size of particles confined in the outer regions is of the order of a centimeter.

## 4  Discussion and conclusion

The proposed model can account for the main characteristics of KH 15D’s observations. It is able to reproduce the light-curve and permits to fit pretty well the data.

(1) The period and duration of the eclipses may be reproduced by the rotation around the star of a giant swarm of solid particles, trapped in a persistent gaseous vortex.

(2) During totality, the rise up of the flux can result from a lower opacity of the central regions, possibly coming from a particle size segregation (the larger the particle, the deeper and the faster the confinement inside the vortex).

(3) The estimated size of the trapped particles ranges from 1 to 10 centimeters, consistent with the grey absorption observed in KH 15D. The presence and growth of bigger particles is very likely due to the high densities and the low relative velocities. This is in agreement with the expected predominance of large grains in PMS disks which seems also required to explain sub-millimeter and millimeter disk observations (eg. Natta et al., 1999).

(4) The shift of the reverse peak observed in successive transits could be explained by gas rotation inside the vortex, the inner zone rotating at half the orbital frequency similarly to what happens in smaller vortices (Chavanis, 2000).

We can also speculate that: (i) $`\mathrm{`}\mathrm{`}`$secular” changes of the light-curve could result from a dynamical evolution of the swarm of particles; (ii) the wavy behaviour at egress could correspond to a turbulent wake of the vortex.

One of the main assumption in our model is that the disk is seen nearly edge-on, at an angle $`i23^{}`$, and with an optical depth $`A_00.60.7`$ magnitude during the out-of-eclipse periods. This is a possible situation if the dust particles have grown and settled down to the mid-plane forming a $`\mathrm{𝑓𝑙𝑎𝑡}`$ sub-disk dominated by large grains. The predominance of large particles in KH 15D is, indeed, completely consistent with our results. However, we want to stress that no numerical simulation is, nowadays, able to reproduce the huge azimutal extent of the vortex required to fit the observations. Further hydrodynamical simulations and theoretical modeling of the particle dynamics are necessary to confirm and improve our model.

Finally, KH 15D appears a good target to test planetary formation models. This is presently the case with the scenario in which vortices can persist in a gas disk acting as traps for dust (protected against losses into the star) and as wombs for planetesimals. Of course, further observations of this object are required to constrain more firmly the model. The possibility to observe large scale structures in circumstellar disks has been recently investigated by Wolf and Klahr (2002), extrapolating the capabilities of the future giant interferometer ALMA. KH 15D and other similar objects are good opportunities for the future space missions like CoRoT, Kepler and Eddington.

A further piece of information was recently provided by Hamilton et al. (2003) using high resolution UVES spectra which confirmed that KH 15D is a weak-lined T Tauri star surrounded by an accretion disk with, possibly, a collimated bipolar jet. Such configuration could favor the accretion-ejection model of Tagger and Pellat (1998) in which persistent vortices can form in the inner region of magnetized disks.

While submitting this letter we learned that H. Klahr also claims the existence of a vortex in the disk of KH 15D to explain the peculiar light-curve of this object. His work seems complementary to ours since he is focusing on the way a giant vortex forms and persists in hydrodynamical simulations, and not on the trapping of solid particles in gaseous structures.

We wish to thank W. Herbst and C.M. Hamilton for communicating us the KH 15D data plotted in Fig. 2. We are also grateful to A. Llebaria for his help in handling the data.


text:
# Comment on “First-principles treatments of electron transport properties for nanoscale junctions”

The recent paper by Fujimoto and Hirose hirose  makes an unfortunate error in discussing the use of the jellium model for the electrodes, which has the effect of making it appear that this model is not adequate to treat the problem of the conductance of gold nanowires. In fact it is entirely adequate, and gives results quite similar to those found in the authors’ more elaborate treatment.

The main point is that the quantity $`D`$ discussed in Sec. III of Ref. hirose, , which represents the distance between the jellium surface (positive-background edge) and the plane of gold atoms (called the “square basis of the nanowires” in Ref. hirose, ) contacting the nanowire is not, as Fujimoto and Hirose suggest, arbitrary or unknown, but by construction of the jellium model, has a perfectly definite value. In their case, this value is (1/4)a0 = 1.93 a.u. (a0 is the gold lattice constant). If the correct value of $`D=1.93`$ a.u. is used in Fig. 6 of Ref. hirose, , which gives conductance vs. $`D`$, then a conductance value of $``$ 0.98 G0 is found, which is to be compared with what the authors call the “true” value of $``$ 1 G0. This shows also that the effect on transport of the interface between the jellium and the gold layer that seems to concern the authors is unimportant.

The essence of the construction of the jellium model is that each lattice plane of the ions of the metal electrode is smeared out symmetrically into a uniform slab of positive charge (see Fig. 1). Thus the plane of ions is at the center of the slab which replaces it,

and the edge of the positive background is half an interplanar spacing in front of the outermost lattice plane (see Fig. 1).

In thinking about this construction, it is useful to consider the difference $`\delta V(𝐫)`$ between the potential due to the ionic lattice and the semi-infinite positive background which represents it. This $`\delta V(𝐫)`$, which can be used to perturbatively re-introduce the discrete lattice as done e.g. in Ref. lang1, , can only be viewed as a small perturbation if the sheet of ions is at the slab center, since otherwise $`\delta V(𝐫)`$ will not vanish at $`\pm \mathrm{}`$ (see footnote 25 of Ref. lang2, ). These issues are also discussed in Ref. lang3, .

Since for Au(100), discussed in Ref. hirose, , the interplanar spacing is (1/2)a0, the distance between the jellium surface and a gold layer put down on it will be half of this, i.e. (1/4)a0. This spacing, as noted before, gives the correct conductance. Similar conclusions can be drawn for other systems as well. For instance, it was shown in Ref. diventra,  that the conductance of an organic molecule between two jellium electrodes is quite similar whether the molecule makes contact directly with the jellium surface or with a plane of gold atoms adjacent to the jellium surface.

One of us (MD) acknowledges support from the NSF Grant Nos. DMR-01-02277 and DMR-01-33075, and Carilion Biomedical Institute.


text:
# An Anisotropic Ballistic Deposition Model with Links to the Ulam Problem and the Tracy-Widom Distribution

## Abstract

We compute exactly the asymptotic distribution of scaled height in a (1+1)–dimensional anisotropic ballistic deposition model by mapping it to the Ulam problem of finding the longest nondecreasing subsequence in a random sequence of integers. Using the known results for the Ulam problem, we show that the scaled height in our model has the Tracy-Widom distribution appearing in the theory of random matrices near the edges of the spectrum. Our result supports the hypothesis that various growth models in $`(1+1)`$ dimensions that belong to the Kardar-Parisi-Zhang universality class perhaps all share the same universal Tracy-Widom distribution for the suitably scaled height variables.

PACS numbers: 05.40.-a, 68.35.Ct, 02.50.-r, 81.10.Aj

Growth processes are ubiquitous in nature. The past few decades have seen an extensive research on a wide variety of both discrete and contiuous growth models Meakin ; KS ; HZ . A large class of these growth models such as the Eden model Eden , restricted solid on solid (RSOS) models RSOS , directed polymers HZ , polynuclear growth models (PNG) PNG  and ballistic deposition models (BD) BD  are believed to belong to the same universality class as that of the Kardar-Parisi-Zhang (KPZ) equation describing the growth of interface fluctuations KPZ . This universality is, however, somewhat restricted in the sense that it refers only to the width or the second moment of the height fluctuations characterized by two independent exponents (the growth exponent $`\beta `$ and the dynamical exponent $`z`$) and the associated scaling function. Moreover, even this restricted universality is established mostly numerically. Only in very few special discrete (1+1)-D models, the exponents $`\beta =1/3`$ and $`z=3/2`$ can be computed exactly via the Bethe ansatz technique Bethe . A natural and important question is whether this universality can be extended beyond the second moment of height fluctuations. For example, does the full distribution of the height fluctuations (suitably scaled) is universal, i.e. is the same for different growth models belonging to the KPZ class? Moreover, the KPZ-type equations are usually attributed to models with small gradients in the height profile and the question whether the models with large gradients belong to the KPZ universality class is still open.

Recently Prähofer and Spohn PS  found an exact mapping between a specific PNG model and the so-called longest increasing subsequence (LIS) problem, also known as the Ulam problem. The LIS problem was first raised by Ulam in the early 60’s Ulam , then the interest in it reappeared in the mathematical literature in 70’s since the work of Vershik and Kerov VK . The exact mapping of PNG to LIS and the subsequent utilization of the exact results available for the LIS problem allowed Prähofer and Spohn to find (besides the exact KPZ growth exponent $`\beta =1/3`$) the exact asymptotic height distribution in the PNG model PS . This distribution turned to be the well known Tracy-Widom distribution appearing in the theory of edge states of random matrices TW . Around the same time, Johansson showed rigorously Johansson  that a specific (1+1)-D directed polymer model, believed to be in the KPZ universality class, also has the same Tracy-Widom distribution for the scaled height (energy) variable. Gravner et. al. found the same Tracy-Widom distribution in another class of $`(1+1)`$ dimensional growth models which they called ‘oriented digital boiling’ modelGTW . It would be interesting to know whether there are other growth models such as the RSOS or the BD ones, which are believed to be in the KPZ universality class as far as the second moment is concerned, also share the same Tracy-Widom distribution for the scaled height.

The purpose of this Letter is to present a BD model which can be mapped exactly to the LIS problem and hence it shares the same Tracy-Widom distribution as the PNG model. This exact result, in combination with the results of PS ; Johansson ; GTW , then lends support to the hypothesis that perhaps all these different growth models, at least in (1+1) dimensions, share the same universal Tracy-Widom distribution for scaled height. This hypothesis, if true, puts the universality on a much stronger footing going beyond just the second moment. Incidentally, to our knowledge, our model is the first exact solution for the full asymptotic height distribution of BD type systems.

Before describing our model, it is worth summarizing the main results for the LIS problem that we use later. Take a set of $`n`$ integers $`\{1,2,3,\mathrm{},n\}`$. Consider all the $`n!`$ possible permutations of this sequence. For any given permutation, let us find out all possible increasing subsequences (terms of a subsequence need not necessarily be consecutive elements) and from them find out the longest one. For example, take $`n=10`$ and consider a particular permutation $`\{8,2,7,\underset{¯}{1},\underset{¯}{3},\underset{¯}{4},10,\underset{¯}{6},\underset{¯}{9},5\}`$. From this sequence, one can form several increasing subsequences such as $`\{8,10\}`$, $`\{2,3,4,10\}`$, $`\{1,3,4,10\}`$ etc. The longest one of all such subsequences is either $`\{1,3,4,6,9\}`$ as shown by the underscores or $`\{2,3,4,6,9\}`$. The length $`l_n`$ of the LIS (in our example $`l_n=5`$) is a random variable as it varies from one permutation to another. In the Ulam problem one considers all the $`n!`$ permutations to be equally likely. Given this uniform measure over the space of permutations, what is the statistics of the random variable $`l_n`$? Ulam found numerically for the average length $`l_n`$ the asymptotic behavior $`l_nc\sqrt{n}`$ for large $`n`$. Later this result was established rigorously by Hammersley Hammersley  and the constant $`c=2`$ was found by Vershik and Kerov VK . Recently, in a seminal paper, Baik, Deift and Johansson (BDJ) BDJ  derived the full distribution of $`l_n`$ for large $`n`$. In particular, they showed that $`l_n2\sqrt{n}+n^{1/6}\chi `$ for large $`n`$, where the random variable $`\chi `$ has an $`n`$-independent distribution which happens to be the Tracy-Widom distribution for the largest eigenvalue of a random matrix drawn from the Gaussian Unitary Ensemble TW . They also showed that when the sequence length $`n`$ itself is a random variable drawn from a Poisson distribution with mean $`n=\lambda `$, the length of the LIS converges for large $`\lambda `$ to

$$l_\lambda 2\sqrt{\lambda }+\lambda ^{1/6}\chi ,$$
(1)
where $`\chi `$ has the Tracy-Widom distribution. The fixed $`n`$ and the fixed $`\lambda `$ ensembles are like the canonical and the grand canonical ensembles in statistical mechanics. The detailed form of the Tracy-Widom distribution is rather complicated and not very illuminating (see PS  for a picture). The BDJ results led to an avalanche of subsequent mathematical works AD .

In our (1+1)-D BD model columnar growth occurs sequentially on a linear substrate consisting of $`L`$ columns with free boundary conditions. The time $`t`$ is discrete and is increased by $`1`$ with every deposition event. One starts at $`t=0`$ with an empty substrate. At any stage of the growth, a column (say the $`k`$-th column) is chosen at random with probability $`p=\frac{1}{L}`$ and a ”brick” is deposited there which increases the height of this column by one unit, $`H_kH_k+1`$. Once this ”brick” is deposited, it screens all the sites at the same level in all the columns to its right from future deposition, i.e. the heights at all the columns to the right of the $`k`$-th column must be strictly greater than or equal to $`H_k+1`$ at all subsequent times. For example, in Fig.1, the first brick (denoted by 1) gets deposited at $`t=1`$ in the 4-th column and it immediately screens all the sites to its right. Then the second brick (denoted by 2) gets deposited at $`t=2`$ again in the same 4-th column whose height now becomes 2 and thus the heights of all the columns to the right of the 4-th column must be $`2`$ at all subsequent times and so on. Formally such growth is implemented by the update rule,

$$H_k(t+1)=\mathrm{max}\{H_k(t),H_{k1}(t),\mathrm{},H_1(t)\}+1.$$
(2)
The model is anisotropic and evidently even the average height profile $`H_k(t)`$ depends nontrivially on both the column number $`k`$ and time $`t`$. Our goal is to compute the asymptotic height distribution $`P_k(H,t)`$ for large $`t`$.

It is easy to find the height distribution $`P_1(H,t)`$ of the first column, since the height there does not depend on any other column. At any stage, the height in the first column either increases by one unit with probability $`p=\frac{1}{L}`$ (if this column is selected for deposit) or stays the same with probability $`1p`$. Thus $`P_1(H,t)`$ is simply the binomial distribution, $`P_1(H,t)=\left(\genfrac{}{}{0pt}{}{t}{H}\right)p^h(1p)^{tH}`$ with $`Ht`$. The average height of the first column thus increases as $`H_1(t)=pt`$ for all $`t`$ and its variance is given by $`\sigma _1^2(t)=tp(1p)`$. While the first column is thus trivial, the dynamics of heights in other columns is nontrivial due to the right-handed infinite range interactions between the columns. For convenience, we subsequently measure the height of any other column with respect to the first one. Namely, by height $`h_k(t)`$ we mean the height difference between the $`(k+1)`$-th column and the first one, $`h_k(t)=H_{k+1}(t)H_1(t)`$, so that $`h_0(t)=0`$ for all $`t`$.

To make progress for columns $`k>0`$, we first consider a (2+1)-D construction of the heap as shown in Fig.2, by adding an extra dimension indicating the time $`t`$. In Fig.2, the $`x`$ axis denotes the column number, the $`y`$ axis stands for the time $`t`$ and the $`z`$ axis is the height $`h`$. In this figure, every time a new block is added, it ”wets” all the sites at the same level to its ”east” (along the $`x`$ axis) and to its ”north” (along the time axis). Here ”wetting” means ”screening” from further deposition at those sites at the same level. This $`(2+1)`$-D system of ”terraces” is in one-to-one correspondence with the $`(1+1)`$-D heap in Fig.1. This construction is reminiscent of the 3D anisotropic directed percolation (ADP) problem studied by Rajesh and Dhar RD . Note however, that unlike the ADP problem, in our case each row labelled by $`t`$ can contain only one deposition eventADP .

The next step is to consider the projection onto the 2D $`(x,y)`$-plane of the level lines separating the adjacent terraces whose heights differ by $`1`$. In this projection, some of the level lines may overlap partially on the plane. To avoid the overlap for better visual purposes, we make a shift $`(x,y)(x+h(x,y),y)`$ and represent these shifted directed lines on the 2D plane in Fig.3.

The black dots in Fig.3 denote the points where the deposition events took place and the integer next to a dot denotes the time of this event. Note that each row in Fig.3 contains a single black dot, i.e. only one deposition per unit of time can occur. In Fig.3, there are 8 such events whose deposition times form the sequence $`\{1,2,3,4,5,6,7,8\}`$ of length $`n=8`$. Now let us read the deposition times of the dots sequentially, but now column by column and vertically from top to bottom in each column, starting from the leftmost one. Then this sequence reads $`\{8,3,5,1,2,6,4,7\}`$ which is just a permutation of the original sequence $`\{1,2,3,4,5,6,7,8\}`$. In the permuted sequence $`\{8,3,5,1,2,6,4,7\}`$ there are $`3`$ LIS’s: $`\{3,5,6,7\}`$, $`\{1,2,6,7\}`$ and $`\{1,2,4,7\}`$, all of the same length $`l_n=4`$. There is a greedy algorithm called the “patience sorting” game devised by Aldous and Diaconis to determine this length of the LISAD . This game goes as follows: start forming piles with the numbers in the permuted sequence starting with the first element which is $`8`$ in our example. So, the number 8 forms the base of the first pile. The next element, if less than 8, goes on top of 8. If not, it forms the base of a new pile. One follows a greedy algorithm: for any new element of the sequence, check all the top numbers on the existing piles starting from the first pile and if the new number is less than the top number of an already existing pile, it goes on top of that pile. If the new number is larger than all the top numbers of the existing piles, this new number forms the base of a new pile. Thus in our example, we form $`4`$ distinct piles: $`[\{8,3,1\},\{5,2\},\{6,4\},\{7\}]`$. The number of piles ($`4`$) is same as the length $`l_n=4`$ of the LIS of this permuted sequence. In fact, Aldous and Diaconis proved AD  that the length of the LIS $`l_n`$ is exactly equal to the number of piles in the corresponding patience sorting game.

Let us note one immediate fact from Fig.3. The numbers belonging to the different level lines in Fig.3 are in one-to-one correspondence with the piles $`[\{8,3,1\},\{5,2\},\{6,4\},\{7\}]`$ in Aldous–Diaconis patience sorting game. Hence, each pile can be identified with an unique level line. Now, the height $`h(x,t)`$ at any given point $`(x,t)`$ in Fig.3 is equal to the number of level lines inside the rectangle bounded by the corners: $`[0,0],[x,0],[0,t],[x,t]`$. Thus, we have the correspondonce: height $``$ number of level lines $``$ number of piles $``$ length $`l_n`$ of the LIS. However, to compute $`l_n`$, we need to know $`n`$ which is the number of black dots inside this rectangle.

Once the problem is reduced to finding the number of black dots or deposition events, we no longer need the Fig.3 (as it may confuse due to the visual shift $`(x,y)(x+h(x,y),y)`$) and can go back to Fig.2, where the north-to-east corners play the same role as the black dots in Fig.2. In Fig.2, to determine the height $`h_k(t)`$ of the $`k`$-th column at time $`t`$, we need to know the number of deposition events inside the $`2`$D plane rectangle $`R_{k,t}`$ bounded by the four corners $`[0,0],[k,0],[0,t],[k,t]`$. Let us begin with the last column $`k=L`$. For $`k=L`$ the number of deposition events $`n`$ in the rectangle $`R_{L,t}`$ is equal to the time $`t`$ because there is only one deposition event per time. In our example $`n=t=8`$. For a general $`k<L`$ the number of deposition events $`n`$ inside the rectangle $`R_{k,t}`$ is a random variable, since some of the rows inside the rectangle may not contain a north-to-east corner or a deposition event. The probability distribution $`P_{k,t}(n)`$ (for a given $`[k,t]`$) of this random variable can, however, be easily found as follows. At each step of deposition, a column is chosen at random from any of the $`L`$ columns. Thus, the probability that a north-to-east corner will fall on the segment of line $`[0,k]`$ (where $`kL`$) is equal to $`k/L`$. The deposition events are completely independent of each other, indicating the absence of correlations between different rows labelled by $`t`$ in Fig.2. So, we are asking the question: given $`t`$ rows, what is the probability that $`n`$ of them will contain a north-to-east corner? This is simply given by the binomial distribution

$$P_{k,t}(n)=\left(\genfrac{}{}{0pt}{}{t}{n}\right)\left(\frac{k}{L}\right)^n\left(1\frac{k}{L}\right)^{tn},$$
(3)
where $`nt`$. Now we are reduced to the following problem: given a sequence of integers of length $`n`$ (where $`n`$ itself is random and is taken from the distribution in Eq.(3)), what is the length of the LIS? Recall that this length is precisely the height $`h_k(t)`$ of the $`k`$-th column at time $`t`$ in our model. In the thermodynamic limit $`L\mathrm{}`$ for $`t1`$ and any fixed $`k`$ such that the quotient $`\lambda =\frac{tk}{L}`$ remains fixed but is arbitrary, the distribution in Eq.(3) becomes a Poisson distribution $`P(n)e^\lambda \frac{\lambda ^n}{n!}`$, with the mean $`\lambda =\frac{tk}{L}`$. We can then directly use the BDJ result in Eq.(1) to predict our main result for the height in the BD model,

$$h_k(t)2\sqrt{\frac{tk}{L}}+\left(\frac{tk}{L}\right)^{1/6}\chi ,$$
(4)
for large $`\lambda =tk/L`$, where the random variable $`\chi `$ has the Tracy-Widom distribution. Using the known exact value $`\chi =1.7711\mathrm{}`$ from the Tracy-Widom distributionTW , we find exactly the asymptotic average height profile in the BD model,

$$h_k(t)2\sqrt{\frac{tk}{L}}1.7721\mathrm{}\left(\frac{tk}{L}\right)^{1/6}.$$
(5)
The leading square root dependence of the profile on the column number $`k`$ has been seen numerically details . The Eq.(5) also predicts an exact sub-leading term with $`k^{1/6}`$ dependence. Similarly, for the variance, $`\sigma _k^2(t)=[h_k(t)h_k(t)]^2`$, we find asymptotically: $`\sigma _k^2(t)c_0\left(\frac{tk}{L}\right)^{1/3}`$, where $`c_0=[\chi \chi ]^2=0.8132\mathrm{}`$ TW . Eliminating the $`t`$ dependence for large $`t`$ between the average and the variance, we get, $`\sigma _k^2(t)ah_k(t)^{2\beta }`$ where the constant $`a=c_0/2^{2/3}=0.51228\mathrm{}`$ and $`\beta =1/3`$, thus recovering the KPZ scaling exponent. In addition to the BD model with infinite range right-handed interaction reported here, we have also analyzed the model (analytically within a mean field theory and numerically) when the right-handed interaction is short rangeddetails . Surprisingly, we found that the asymptotic average height profile is independent of the range of interactiondetails .

In summary, we have shown that the asymptotic scaled height in an anisotropic (1+1)D BD model has the Tracy-Widom distribution. Our exact result, in combination with those of Refs.PS ; Johansson ; GTW  where the same distribution was found in rather different growth models, suggests that the universality in all these growth processes is perhaps much wider than it was thought before, extending to the full asymptotic height distribution and is not restricted only to the second moment of the height.

We thank the hospitality of the Institut Henri Poincaré in Paris where this work was initiated during the trimester ”Geometry and Statistics of Random Growth” in 2003.


text:
# Transverse momentum distribution of net baryon number at RHIC

## References


text:
# Time evolution, cyclic solutions and geometric phases for general spin in an arbitrarily varying magnetic field

thanks: published in J. Phys. A 36 (2003) 6799-6806.thanks: not for correspondence

## Abstract

A neutral particle with general spin and magnetic moment moving in an arbitrarily varying magnetic field is studied. The time evolution operator for the Schrödinger equation can be obtained if one can find a unit vector that satisfies the equation obeyed by the mean of the spin operator. There exist at least $`2s+1`$ cyclic solutions in any time interval. Some particular time interval may exist in which all solutions are cyclic. The nonadiabatic geometric phase for cyclic solutions generally contains extra terms in addition to the familiar one that is proportional to the solid angle subtended by the closed trace of the spin vector.

Since the discovery of the geometric phase berry ; simon ; aha ; sam ; wu-li ; jordan ; jackiw ; resource ; li-book , particles with spin and magnetic moment moving in time-dependent magnetic fields have received much attention wang ; wagh1 ; wagh2 ; fer ; fer-pla ; layton ; gao ; ni ; zhang ; zhu00 ; pra01 ; jpa01 ; jpa02 ; ni-book , though the subject is rather old and some discussions can be found in the textbook landau . Neutral particles are of special interest since the problem is easier and the Schrödinger equation can be solved analytically in some special cases, say, uniform magnetic fields with a fixed direction or rotating ones. Thus the model is very suitable for the study of time evolution, cyclic solutions and geometric phases etc. However, some problems in this model are still not clear. First, the Schrödinger equation for general spin, or even for spin $`1/2`$, in an arbitrarily varying magnetic field seems impossible to be solved analytically. Second, though the existence of cyclic solutions in a given time interval may be ensured by the existence of eigenvectors for the unitary time evolution operator, it seems not clear there are how many ones in the general case. Third, for spin $`1/2`$ it is well known that any cyclic solution in an arbitrary magnetic field has a nonadiabatic geometric phase proportional to the solid angle subtended by the closed trace of the spin vector. For higher spin, however, this is true only for cyclic solutions with special initial conditions layton ; gao ; jpa02 . For more general cyclic solutions in a rotating magnetic field, we have shown that the nonadiabatic geometric phase contains an extra term in addition to the one proportional to the solid angle. The extra term vanishes automatically for spin $`1/2`$. For higher spin, however, it depends on the initial condition jpa02 . It is still not clear what is the relation between the nonadiabatic geometric phase and the solid angle for general cyclic solutions in an arbitrary magnetic field. In this paper we are going to deal with these problems, and try to solve them to some extent. Besides the theoretical interest in itself and other applications li-book ; zhu00 , this subject has been recently recognized to be of great interest in the physics of quantum computation wang01 ; zhu02 .

Consider a neutral particle with spin $`s`$ ($`s=1/2,1,3/2,\mathrm{}`$) and magnetic moment $`𝝁=\mu 𝐬/s`$, where $`𝐬`$ is the spin operator in the unit of $`\mathrm{}`$, satisfying $`[s_i,s_j]=\mathrm{i}ϵ_{ijk}s_k`$ (for spin $`1/2`$, $`𝐬=𝝈/2`$ and $`𝝁=\mu 𝝈`$). In a uniform but time-dependent magnetic field $`𝐁(t)=B(t)𝐧(t)`$ where $`𝐧(t)`$ is a unit vector, it has the Hamiltonian $`H(t)=𝝁𝐁(t)=\mathrm{}\omega _B(t)𝐬𝐧(t)`$, where $`\omega _B(t)=\mu B(t)/s\mathrm{}`$, and the Schrödinger equation $`\mathrm{i}\mathrm{}_t\psi (t)=H(t)\psi (t)`$ takes the form

$$_t\psi (t)=\mathrm{i}\omega _B(t)𝐬𝐧(t)\psi (t).$$
(1)
Define the spin vector as

$$𝐯(t)=\mathbf{(}\psi (t),𝐬\psi (t)\mathbf{)}.$$
(2)
Using Eq. (1) it is easy to show that it obeys the equation

$$\dot{𝐯}(t)=\omega _B(t)𝐧(t)\times 𝐯(t).$$
(3)
We are not going to solve Eq. (1) in the general case since this seems impossible. However, we will show that the time evolution operator for Eq. (1) can be obtained without any chronological product if one can find one nontrivial (nonzero) solution, say, a unit vector $`𝐞(t)`$, to Eq. (3). This is of interest since the latter is easier and more cases can be solved fer-pla . Actually, Eq. (1) involves operators while Eq. (3) involves only c-numbers. On the other hand, if transformed to a matrix equation, Eq. (1) involves $`2s+1`$ complex variables, while Eq. (3) involves only three real ones (actually two since it is easy to see that $`𝐯^2(t)=𝐯^2(0)`$). Using the time evolution operator and Eq. (3), one can discuss cyclic solutions and geometric phases in a most general way. In particular, we will show that there exist at least $`2s+1`$ cyclic solutions in any time interval. A general relation between the nonadiabatic geometric phase and the solid angle subtended by the closed trace of the spin vector will be established.

To begin, we take an arbitrary unit vector $`𝐞_0`$, and the eigenstate of $`𝐬𝐞_0`$ with eigenvalue $`m`$ will be denoted by $`\chi _m`$. We take the initial state of the system to be $`\psi (0)=\chi _m`$, that is

$$𝐬𝐞_0\psi (0)=m\psi (0),m=s,s1,\mathrm{},s.$$
(4)
Obviously $`𝐯(0)=m𝐞_0`$ in this initial state. Now we define a vector $`𝐞(t)`$ by Eq. (3) with the initial condition $`𝐞_0`$, that is

$$\dot{𝐞}(t)=\omega _B(t)𝐧(t)\times 𝐞(t)$$
(5)
with $`𝐞(0)=𝐞_0`$. We would assume that $`𝐁(t)`$ varies continuously, so that any solution $`𝐞(t)`$ is well behaved. As pointed out above, $`𝐞^2(t)=𝐞_0^2`$, so $`𝐞(t)`$ is a unit vector at any time. We have proven in Ref. jpa02  that

$$𝐬𝐞(t)\psi (t)=m\psi (t)$$
(6)
holds at all later times. To be self-contained, we repeat here the proof by induction.

By definition, Eq. (6) is valid at $`t=0`$. We assume that it is valid at time $`t`$, what we need to do is to show that it is also true at time $`t+\mathrm{\Delta }t`$ where $`\mathrm{\Delta }t`$ is an infinitesimal increment of time. In fact, using Eqs. (1) and (5) we have

$$\psi (t+\mathrm{\Delta }t)=\psi (t)+\mathrm{i}\omega _B(t)𝐬𝐧(t)\psi (t)\mathrm{\Delta }t,$$
(7a)
$$𝐞(t+\mathrm{\Delta }t)=𝐞(t)\omega _B(t)𝐧(t)\times 𝐞(t)\mathrm{\Delta }t.$$
(7b)
After some simple algebra, the conclusion is achieved.

Because $`𝐞(t)`$ is a unit vector, we can write in some rectangular coordinates

$$𝐞(t)=\mathbf{(}\mathrm{sin}\theta (t)\mathrm{cos}\varphi (t),\mathrm{sin}\theta (t)\mathrm{sin}\varphi (t),\mathrm{cos}\theta (t)\mathbf{)}.$$
(8)
Using the formula jpa02 ; wilcox

$$\mathrm{e}^{\mathrm{i}\xi 𝐬𝐛}𝐬\mathrm{e}^{\mathrm{i}\xi 𝐬𝐛}=[𝐬(𝐬𝐛)𝐛]\mathrm{cos}\xi +(𝐛\times 𝐬)\mathrm{sin}\xi +(𝐬𝐛)𝐛,$$
(9)
where $`𝐛`$ is any unit vector, it is not difficult to show that

$$𝐬𝐞(t)=\mathrm{e}^{\mathrm{i}\theta (t)𝐬𝐝(t)}s_z\mathrm{e}^{\mathrm{i}\theta (t)𝐬𝐝(t)},$$
(10)
where

$$𝐝(t)=\mathbf{(}\mathrm{sin}\varphi (t),\mathrm{cos}\varphi (t),0\mathbf{)}.$$
(11)
Therefore the eigenstate of $`𝐬𝐞(t)`$ with eigenvalue $`m`$ is

$$\psi (t)=\mathrm{e}^{\mathrm{i}\alpha _m(t)}\mathrm{e}^{\mathrm{i}\theta (t)𝐬𝐝(t)}\chi _m^0,$$
(12)
where $`\chi _m^0`$ is the eigenstate of $`s_z`$ with eigenvalue $`m`$, and $`\alpha _m(t)`$ is a phase that cannot be determined by the eigenvalue equation. However, $`\alpha _m(t)`$ is not arbitrary. To satisfy the Schrödinger equation, it should be determined by the other variables $`\theta (t)`$ and $`\varphi (t)`$. In fact, the above equation yields

$$\mathbf{(}\psi (t),\psi (t+\mathrm{\Delta }t)\mathbf{)}=1+\mathrm{i}\dot{\alpha }_m(t)\mathrm{\Delta }t+\mathbf{(}\chi _m^0,\mathrm{e}^{\mathrm{i}\theta (t)𝐬𝐝(t)}_t\mathrm{e}^{\mathrm{i}\theta (t)𝐬𝐝(t)}\chi _m^0\mathbf{)}\mathrm{\Delta }t.$$
(13)
Using the formula wilcox

$$\mathrm{e}^{F(t)}_t\mathrm{e}^{F(t)}=_0^1\mathrm{e}^{\lambda F(t)}\dot{F}(t)e^{\lambda F(t)}d\lambda ,$$
(14)
where $`F(t)`$ is any operator depending on $`t`$, and then using Eq. (9), we obtain

$$\mathbf{(}\psi (t),\psi (t+\mathrm{\Delta }t)\mathbf{)}=1+\mathrm{i}\dot{\alpha }_m(t)\mathrm{\Delta }t+\mathrm{i}m[1\mathrm{cos}\theta (t)]\dot{\varphi }(t)\mathrm{\Delta }t.$$
(15)
On the other hand, from Eq. (1) we have

$$\mathbf{(}\psi (t),\psi (t+\mathrm{\Delta }t)\mathbf{)}=1+\mathrm{i}\omega _B(t)𝐯(t)𝐧(t)\mathrm{\Delta }t.$$
(16)
Note that $`𝐯(t)`$ and $`𝐞(t)`$ satisfy the same equation, and $`𝐯(0)=m𝐞_0`$, we have $`𝐯(t)=m𝐞(t)`$. Comparing the two results above and taking this relation into account, we obtain

$$\dot{\alpha }_m(t)=m[1\mathrm{cos}\theta (t)]\dot{\varphi }(t)+m\omega _B(t)𝐞(t)𝐧(t).$$
(17)
Therefore

$$\alpha _m(t)\alpha _m(0)=m\alpha (t),$$
(18)
where

$$\alpha (t)=_0^t[1\mathrm{cos}\theta (t^{})]\dot{\varphi }(t^{})dt^{}+_0^t\omega _B(t^{})𝐞(t^{})𝐧(t^{})dt^{}.$$
(19)
Substituting into Eq. (12) we obtain

$$\psi (t)=\mathrm{e}^{\mathrm{i}\theta (t)𝐬𝐝(t)}\mathrm{e}^{\mathrm{i}\alpha (t)s_z}\mathrm{e}^{\mathrm{i}\theta (0)𝐬𝐝(0)}\psi (0).$$
(20)
We denote the time evolution operator as $`U(t)`$, defined by the equation $`\psi (t)=U(t)\psi (0)`$ with an arbitrary $`\psi (0)`$, then the above equation is equivalent to

$$U(t)\chi _m=\mathrm{e}^{\mathrm{i}\theta (t)𝐬𝐝(t)}\mathrm{e}^{\mathrm{i}\alpha (t)s_z}\mathrm{e}^{\mathrm{i}\theta (0)𝐬𝐝(0)}\chi _m.$$
(21)
Now an arbitrary initial state $`\psi (0)`$ can always be expanded as

$$\psi (0)=\underset{m}{}c_m\chi _m.$$
(22)
Applying $`U(t)`$ to both sides of this equation, using Eq. (21), and noting that the operators on the right-hand side of that equation is independent of $`m`$, we immediately realize that Eq. (20) is in fact valid for an arbitrary initial state. Thus we arrive at the result

$$U(t)=\mathrm{e}^{\mathrm{i}\theta (t)𝐬𝐝(t)}\mathrm{e}^{\mathrm{i}\alpha (t)s_z}\mathrm{e}^{\mathrm{i}\theta (0)𝐬𝐝(0)}.$$
(23a)
Using Eq. (9), it can be recast in the form
$$U(t)=\mathrm{e}^{\mathrm{i}\theta (t)𝐬𝐝(t)}\mathrm{e}^{\mathrm{i}\theta (0)𝐬𝐝(0)}\mathrm{e}^{\mathrm{i}\alpha (t)𝐬𝐞_0}.$$
(23b)
Eq. (23b) is suitable for the general discussions below while Eq. (23a) may be more convenient for practical calculations.

Let us make some remarks on the result. First, we see that once a solution of Eq. (5) is found, the time evolution operator for Eq. (1) is available and it involves no chronological product. The result depends formally on $`𝐞_0`$, but $`𝐞_0`$ is merely an auxiliary object, hence the result must be essentially independent of it, though it might be difficult to prove this explicitly. In practical calculations, one should choose a solution $`𝐞(t)`$ that is as simple as possible such that $`U(t)`$ can be easily reduced to the simplest form. When this approach is used to the simple cases such as rotating magnetic fields or ones with a fixed direction, it indeed leads to the same results as those obtained previously. Second, the operator $`U(t)`$ depends not only on $`𝐞(t)`$, but also on the history of it. This is obvious from Eq. (19). Third, though $`\varphi (t)`$ is indefinite when $`\theta (t)=0`$ or $`\pi `$, the above result $`U(t)`$ is in fact well behaved everywhere. That it is well defined at $`\theta (t)=0`$ is obvious. If $`\theta (0)=0`$, there is no problem either. The case with $`\theta (0)=\pi `$ can be avoided since one can always choose a coordinate system such that $`\theta (0)\pi `$. However, for a general evolution, the case with $`\theta (t)=\pi `$ at some instant cannot be avoided. Thus we must show that $`U(t)`$ is well behaved at $`\theta (t)=\pi `$. Suppose that $`\theta (t_0)=\pi `$, then we have $`\varphi (t_0^+)\varphi (t_0^{})=\pm \pi `$, $`𝐝(t_0^+)=𝐝(t_0^{})`$, and $`\alpha (t_0^+)\alpha (t_0^{})=2\pi `$. With these relations it is not difficult to show that $`U(t_0^+)=U(t_0^{})`$. Since both $`U(t_0^+)`$ and $`U(t_0^{})`$ are well defined, we may define $`U(t_0)=lim_{tt_0}U(t)`$. This makes $`U(t)`$ well defined and continuous at $`t=t_0`$. Fourth, by straightforward calculations it can be shown that $`_tU(t)=\mathrm{i}\omega _B(t)𝐬𝐧(t)U(t)`$ and $`U(0)=1`$, as expected.

Now we can go further to discuss cyclic solutions in any time interval $`[0,\tau ]`$ where $`\tau `$ is an arbitrarily given time. These cyclic solutions are not necessarily cyclic in subsequent time intervals with the same length, say, $`[\tau ,2\tau ]`$.

Since Eq. (5) is a linear differential equation, the general solution $`𝐞(t)`$ must depend on the initial vector $`𝐞_0`$ linearly. Thus it can be written in a matrix form

$$e_i(t)=E_{ij}(t)e_{0j},$$
(24)
where the matrix $`E(t)`$ is obviously real. If both $`𝐞_1(t)`$ and $`𝐞_2(t)`$ are solutions to Eq. (5), it is easy to show that $`𝐞_1(t)𝐞_2(t)=𝐞_1(0)𝐞_2(0)`$. Therefore the matrix $`E(t)`$ is an orthogonal one, and its eigenvalues at any time $`t`$ has the form $`\{1,\sigma (t),\sigma ^{}(t)\}`$, where $`\sigma (t)`$ is a complex number with $`|\sigma (t)|=1`$, and $`\sigma ^{}(t)`$ its complex conjugate.

If $`\sigma (\tau )1`$, one eigenvector $`𝜼(\tau )`$ of the matrix $`E(\tau )`$ with eigenvalue $`1`$ can be found, which satisfies $`E_{ij}(\tau )\eta _j(\tau )=\eta _i(\tau )`$. It can be taken as real and normalized. Now if we choose

$$𝐞_0=𝜼(\tau ),$$
(25)
we have $`e_i(\tau )=E_{ij}(\tau )e_{0j}=E_{ij}(\tau )\eta _j(\tau )=\eta _i(\tau )=e_i(0)`$, that is

$$𝐞(\tau )=𝐞_0.$$
(26)
This means that $`\theta (\tau )=\theta (0)`$ and $`𝐝(\tau )=𝐝(0)`$, and leads to

$$U(\tau )=\mathrm{e}^{\mathrm{i}\alpha (\tau )𝐬𝐞_0}.$$
(27)
Now it is clear that with the initial condition $`\psi (0)=\chi _m`$ ($`m=s,s1,\mathrm{},s`$), we have a cyclic solution in the time interval $`[0,\tau ]`$. More specifically, $`\psi (\tau )=\mathrm{e}^{\mathrm{i}\delta }\psi (0)`$, where the total phase change is $`\delta =m\alpha (\tau ),\mathrm{mod}2\pi `$, with $`\alpha (\tau )`$ given by

$$\alpha (\tau )=\mathrm{\Omega }_e+_0^\tau \omega _B(t)𝐞(t)𝐧(t)dt,$$
(28)
where

$$\mathrm{\Omega }_e=_0^\tau [1\mathrm{cos}\theta (t)]\dot{\varphi }(t)dt$$

is the solid angle subtended by the closed trace of $`𝐞(t)`$. Notice that $`𝐯(t)=m𝐞(t)`$, the dynamic phase $`\beta =\mathrm{}^1_0^\tau H(t)dt`$ turns out to be

$$\beta =m_0^\tau \omega _B(t)𝐞(t)𝐧(t)dt.$$
(29)
Therefore the nonadiabatic geometric phase is

$$\gamma =\delta \beta =m\mathrm{\Omega }_e,\mathrm{mod}2\pi .$$
(30)
Since $`\mathrm{\Omega }_e=ϵ(m)\mathrm{\Omega }_v,\mathrm{mod}4\pi `$, where $`\mathrm{\Omega }_v`$ is the solid angle subtended by the closed trace of the spin vector, we have finally

$$\gamma =|m|\mathrm{\Omega }_v,\mathrm{mod}2\pi ,$$
(31)
in accord with the results previously obtained layton ; gao ; jpa02 . Thus we see that $`2s+1`$ cyclic solutions are available in any time interval $`[0,\tau ]`$, and all phases can be expressed in terms of the unit vector $`𝐞(t)`$.

States with initial condition other than the above ones are in general not cyclic ones, even those in which $`𝐯(0)`$ points in the direction of $`𝐞_0`$ such that $`𝐯(\tau )=𝐯(0)`$. However, if $`\alpha (\tau )/\pi `$ happens to be a rational number other than an even integer (the case with $`\alpha (\tau )/\pi `$ an even integer will be discussed below), some other cyclic solutions may be available. To be more specific, let $`\alpha (\tau )=p\pi /n`$, where $`n`$ is a natural number and $`p`$ an integer. When $`n=1`$, $`p`$ is an odd number, and when $`n>1`$ it is prime with $`p`$. If $`sn`$, we have cyclic solutions with initial condition, say (no such solution exits for $`s=1/2`$),

$$\psi (0)=\underset{j=N_1}{\overset{N_2}{}}c_j\chi _{m+2nj},$$
(32)
where $`N_1`$ and $`N_2`$ are nonnegative integers, and

$$m2nN_1s,m+2nN_2s,\underset{j=N_1}{\overset{N_2}{}}|c_j|^2=1.$$

In this initial state

$$𝐯(0)=v_0𝐞_0,$$
(33)
where $`v_0=_{j=N_1}^{N_2}(m+2nj)|c_j|^2`$ may be either positive or negative. It is easy to show that $`\psi (\tau )=\mathrm{e}^{\mathrm{i}\delta }\psi (0)`$, where $`\delta =m\alpha (\tau ),\mathrm{mod}2\pi `$. Because of Eq. (33), we have $`𝐯(t)=v_0𝐞(t)`$, and the dynamic phase is

$$\beta =\mathrm{}^1_0^\tau H(t)dt=_0^\tau \omega _B(t)𝐯(t)𝐧(t)dt=v_0_0^\tau \omega _B(t)𝐞(t)𝐧(t)dt.$$

Using Eq. (28), we have

$$\beta =v_0[\alpha (\tau )+\mathrm{\Omega }_e].$$
(34)
This holds regardless of the values of $`v_0`$ and $`\alpha (\tau )`$, as long as Eq. (33) is valid. Suppose that in the process from $`t=0`$ to $`t=\tau `$, $`𝐞(t)`$ encircles the polar axis $`K`$ times ($`K>0`$ for anticlockwise traces and $`K<0`$ for clockwise ones), then we have

$$\mathrm{\Omega }_e+\mathrm{\Omega }_e=4\pi K.$$
(35)
This leads to $`v_0\mathrm{\Omega }_e=|v_0|\mathrm{\Omega }_v+2\pi K(v_0|v_0|)`$. The geometric phase turns out to be

$$\gamma =|v_0|\mathrm{\Omega }_v+2\pi K(|v_0|v_0)+(mv_0)p\pi /n,\mathrm{mod}2\pi .$$
(36)
In this case we see that the geometric phase contains extra terms in addition to the one proportional to $`\mathrm{\Omega }_v`$, unless the sum of these extra terms happens to be an integral multiple of $`2\pi `$. Note that the above relation holds for the case $`v_0=0`$ as well, though $`\mathrm{\Omega }_v`$ is not well defined in this case.

If it happens that $`\alpha (\tau )=2k\pi `$, then $`U(\tau )=\mathrm{e}^{\mathrm{i2}\pi ks}`$ becomes a c-number and all solutions are cyclic in the time interval $`[0,\tau ]`$. However, this is true only on the premise of (26). Thus $`\alpha (\tau )=2k\pi `$ alone is not a sufficient condition for all solutions to be cyclic, but it is easy to show that it is a necessary one. A sufficient condition is $`\sigma (\tau )=1`$.

Now if $`\sigma (\tau )=1`$, $`E(\tau )`$ becomes a unit matrix. In this case any vector is its eigenvector with eigenvalue $`1`$. Therefore Eqs. (26) and (27) hold for any unit vector $`𝐞_0`$. In particular, we have

$$U(\tau )=\mathrm{e}^{\mathrm{i}\alpha _1(\tau )s_x}=\mathrm{e}^{\mathrm{i}\alpha _2(\tau )s_y}=\mathrm{e}^{\mathrm{i}\alpha _3(\tau )s_z},$$
(37)
where $`\alpha _1(\tau )`$ is given by Eq. (28) with $`𝐞_0=𝐞_x`$ and similarly for $`\alpha _2(\tau )`$ and $`\alpha _3(\tau )`$. Since $`s_x`$, $`s_y`$ and $`s_z`$ are independent operators, the above equation cannot hold unless $`\alpha _i(\tau )=2\pi k_i`$ where $`k_i`$ are integers such that $`U(\tau )`$ becomes a c-number. In general, with an initial unit vector $`𝐞_0`$, we have

$$\alpha (\tau )=2\pi k,$$
(38)
and

$$U(\tau )=\mathrm{e}^{\mathrm{i}\alpha (\tau )𝐬𝐞_0}=\mathrm{e}^{\mathrm{i2}\pi ks}.$$
(39)
Let us see what is the dependence of $`\alpha (\tau )`$ or $`k`$ on the direction of $`𝐞_0`$. Consider two initial unit vectors $`𝐞_0`$ and $`𝐞_0^{}`$, whose difference $`\delta 𝐞_0=𝐞_0^{}𝐞_0`$ is infinitesimal (then $`\delta 𝐞_0𝐞_0=0`$). The difference in $`\alpha (\tau )`$ is, according to Eq. (28),

$$\delta \alpha (\tau )=\alpha ^{}(\tau )\alpha (\tau )=\delta \mathrm{\Omega }_e+_0^\tau \omega _B(t)\delta 𝐞(t)𝐧(t)dt.$$
(40)
Since $`e_i(t)=E_{ij}(t)e_{0j}`$, we have $`\delta e_i(t)=E_{ij}(t)\delta e_{0j}`$, and the second term in the above equation is consequently infinitesimal. Moreover, the trace of $`𝐞^{}(t)`$ is very close to that of $`𝐞(t)`$, thus the difference in the solid angles subtended by them is infinitesimal. Therefore both $`\delta \alpha (\tau )`$ and $`\delta k=\delta \alpha (\tau )/2\pi `$ are infinitesimal as well. Now that $`k`$ can take only integer values, an obvious consequence is that $`\delta \alpha (\tau )=0`$ and $`k^{}=k`$.

A subtle case has been overlooked in the above discussions, however. This happens when $`𝐞(t)`$ goes by the south pole ($`\theta =\pi `$) on one side and $`𝐞^{}(t)`$ goes by it on the other. In other words, the closed trace of $`𝐞(t)`$ as a whole goes across the south pole when $`𝐞_0`$ varies to $`𝐞_0^{}`$. In this case $`\alpha (\tau )`$ will be changed by an integral multiple of $`4\pi `$ and $`k`$ by an even integer. One can of course removing this change by rotating the coordinate system such that the trace of $`𝐞(t)`$ does not go across the south pole. However, this can only be done locally, and globally it is impossible in general. In other words, one cannot choose a coordinate system such that $`k`$ is the same for all $`𝐞_0`$, except for some simple cases, say, magnetic fields with a fixed direction. The above dependence of $`\alpha (\tau )`$ on $`𝐞_0`$ has no consequence on the time evolution operator as expected. This is easily seen from Eq. (39): when $`k`$ changes by an even integer, $`U(\tau )`$ remains the same.

In the special case where $`k`$ is the same for all $`𝐞_0`$, a geometric explanation of $`k`$ is available. If we take the unit vector $`𝐞_0^{}=𝐞_0`$ as an initial condition to Eq. (5), then the solution is $`𝐞^{}(t)=𝐞(t)`$. Since $`\alpha (\tau )`$ corresponding to $`𝐞(t)`$ and $`\alpha ^{}(\tau )`$ corresponding to $`𝐞^{}(t)`$ are equal as assumed, we have, according to Eq. (28), $`2\alpha (\tau )=(\mathrm{\Omega }_e+\mathrm{\Omega }_e)`$. As before, $`\mathrm{\Omega }_e+\mathrm{\Omega }_e=4\pi K`$, and $`2\alpha (\tau )=4\pi k`$, so that $`k=K`$. Here $`K`$ is the winding number of $`𝐞(t)`$ around the polar axis. Unfortunately, the geometric meaning of $`k`$ in the general case is still not clear.

Now that $`U(\tau )`$ is a c-number, all solutions to Eq. (1) become cyclic in the time interval $`[0,\tau ]`$. The total phase change is $`\delta =2\pi ks=s\alpha (\tau ),\mathrm{mod}2\pi `$. Since this phase is determined only up to an integral multiple of $`2\pi `$, the dependence of $`\alpha (\tau )`$ or $`k`$ on the initial vector $`𝐞_0`$ does not affect the result. For convenience we take the $`\alpha (\tau )`$ or $`k`$ that is the same as that appears below. If in the initial state $`𝐯(0)0`$, we take $`𝐞_0=𝐯(0)/v_0`$ where $`v_0=|𝐯(0)|>0`$. Thus $`𝐯(t)=v_0𝐞(t)`$ and $`\mathrm{\Omega }_e=\mathrm{\Omega }_v`$. The dynamic phase is, as shown in Eq. (34), $`\beta =v_0[\alpha (\tau )+\mathrm{\Omega }_v]`$. It should be remarked that both $`\alpha (\tau )`$ and $`\mathrm{\Omega }_v`$ may be different by an integral multiple of $`4\pi `$ in different coordinate systems. However, $`\beta `$ has a definite value, independent of the choice of coordinate systems, which can be easily seen from the expression $`\beta =v_0_0^\tau \omega _B(t)𝐞(t)𝐧(t)dt`$. The geometric phase turns out to be

$$\gamma =v_0\mathrm{\Omega }_v+(sv_0)2\pi k,\mathrm{mod}2\pi .$$
(41)
Here the first term is the familiar one, but an extra term appears. If $`s=1/2`$, we have $`v_0=1/2`$ for any initial state, and the above result reduces to $`\gamma =\mathrm{\Omega }_v/2`$, a well-known result. For higher spin, however, the extra term depends on the initial state, as both $`v_0`$ and $`k`$ depend on it. It vanishes (mod $`2\pi `$ of course) when $`sv_0`$ is an integer, especially when the initial state is an eigenstate of $`𝐬𝐞_0`$ (it cannot be an eigenstate of $`𝐬𝐞_0^{}`$ with some other unit vector $`𝐞_0^{}`$ since $`𝐯(0)`$ points in the direction of $`𝐞_0`$), as expected. If in the initial state $`𝐯(0)=0`$, we have $`𝐯(t)=0`$ and thus $`\beta =0`$. It is easy to see that Eq. (41) holds in this case as well, though $`\mathrm{\Omega }_v`$ is not well defined.

The general result (41) has been confirmed by practical calculations in the case of a rotating magnetic field, where both $`\gamma `$ and $`\mathrm{\Omega }_v`$ can be calculated explicitly jpa02 . For a magnetic field with a fixed direction, it is easier to carry out similar calculations to verify this result.

In summary, we have shown that the time evolution operator for the Schrödinger equation (1) can be obtained if one nontrivial solution to Eq. (5) can be found. We proved that at least $`2s+1`$ cyclic solutions of the Schrödinger equation exist in any time interval. These cyclic solutions can be worked out in principle if the general solution to Eq. (5) is known. There may exist some particular time interval in which all solutions are cyclic. The nonadiabatic geometric phase for cyclic solutions contains in general extra terms in addition to the familiar one that is proportional to the solid angle subtended by the trace of the spin vector. For spin $`1/2`$ there is no such extra term.

###### Acknowledgements.

This work was supported by the National Natural Science Foundation of China (10275098), and by the Foundation of the Advanced Research Center of Sun Yat-Sen University (02P3).


text:
# Constant mean curvature slices in the extended Schwarzschild solution and collapse of the lapse. Part II

## A Introduction

This is a sequel of our previous work on the constant mean curvature (CMC) slices of the extended Schwarzschild geometry. Here we get CMC foliations by solving Einstein equations in a particular gauge. A crucial role is played by a condition (Eq. (5) below) imposed on the lapse. While this method is completely equivalent to the other more geometric approach (see ), it seems to be more straightforward and technically simpler. We focus on the concise derivation of the explicit CMC foliation near the critical point of the CMC foliation. The final result is identical to the result derived in .

The constant mean curvature foliations have been recently investigated numerically in the simulation of a single spherically-symmetric black hole . We hope that our analytic results appear helpful in the verification of the numerical schemes.

## B CMC slicing of the Schwarzschild spacetime

The notation is the same as in the preceding paper . We define

$$(pR)^2=4\left[1\frac{2m}{R}+\left(\frac{KR}{3}\frac{C}{R^2}\right)^2\right],$$
(1)
$$\gamma (R,t)=1+8_tC_R^{\mathrm{}}𝑑r\frac{1}{r^5p^3}.$$
(2)
and

$$N=\gamma \frac{pR}{2}.$$
(3)
Here $`m`$ is the mass, $`K`$ (the trace of the extrinsic curvature) is a constant and C is a time-dependent parameter which measures the transverse part of the extrinsic curvature.

The Schwarzschild line element, expressed in terms of coordinates adapted to the constant mean curvature foliation, is given by

$`ds^2`$ $`=`$ $`dt^2\left(N^2\gamma ^2\left({\displaystyle \frac{KR}{3}}{\displaystyle \frac{C}{R^2}}\right)^2\right)+4N{\displaystyle \frac{\frac{C}{R^3}\frac{K}{3}}{p^2R}}dtdR+{\displaystyle \frac{4}{(pR)^2}}dR^2+R^2d\mathrm{\Omega }^2.`$ (4)
The hypersurfaces of constant time are CMC slices, asymptotic to the CMC slices of Minkowskian geometry.

## C Elliptic slicing condition

A minimal surface is a locus of points defined by the condition $`p=0`$. Choose a CMC Cauchy hypersurface $`\mathrm{\Sigma }_C`$ of the extended Schwarzschild manifold corresponding to a parameter $`C`$ and let $`R_0`$ be an areal radius corresponding to a simple zero of $`p^2`$; that is $`p^2(R_0)=0`$ but $`_Rp^2|_{R_0}0`$. Futhermore, assume that

$$\frac{_rN}{\sqrt{a}}|_{R_0}=0$$
(5)
at $`R_0`$. The condition (5) yields

$$_tC=\frac{1}{8I(R_0)}.$$
(6)
Here

$$I(R_0)_{R_0}\frac{dr}{pr}\frac{6\frac{C^2}{r^4}+\frac{K^2r^2}{3}}{\left(2m+\frac{2KC}{3}+\frac{2K^2r^2}{9}\frac{4C^2}{r^3}\right)^2}.$$
(7)
The value of the lapse function $`N`$ at the minimal surface, that is at the areal radius $`R_0`$, can be shown to be equal (using Eqs. (13)) to

$$N=\frac{dC}{dt}\frac{1}{m+\frac{KC}{3}+\frac{K^2R_0^3}{9}2\frac{C^2}{R_0^3}}.$$
(8)
The lapse $`N`$ is strictly positive at the minimal surface corresponding to a simple zero $`R_0`$. Eqs. (13) imply that $`N(R)>N(R_0)`$ if $`R>R_0`$ and therefore the lapse exists on all of $`\mathrm{\Sigma }_C`$. Equation (6) dictates the rate of change of the parameter $`C`$. It is clear that one can uniquely construct a foliation of a part of the extended Schwarzschild geometry by imposing the condition (5) at minimal surfaces on all slices to the future of a given one. The leaves of the resulting foliation connect two null infinities of the extended Schwarzschild spacetime. This gives us a curve $`R_0(t)`$ of zeroes of the mean curvature $`p`$. It is evident, just by inspecting the explicit solution presented above, that the line running along the locations of minimal surfaces $`R_0(t)`$ can be arranged to be smooth. It can be chosen to coincide with the ‘vertical’ $`t=0`$ axis in standard Schwarzschild coordinates.

This construction breaks down when $`R_0`$ ceases to be a simple zero of $`p^2`$, since expressions appearing in Eqs. (6) and (8) become unboundedly large. The goal of this paper is to show the asymptotic behaviour of the lapse at the critical minimal surface.

## D The evolution of C near critical point

Let $`C_{}`$ and $`R_{}`$ be degenerate, that is such that the zero of $`p^2`$ ceases to be simple. In this case both $`p`$ and its derivative $`_Rp`$ vanish; that means that

$`1{\displaystyle \frac{2m}{R_{}}}{\displaystyle \frac{2KC_{}}{3R_{}}}+{\displaystyle \frac{K^2R_{}^2}{9}}+{\displaystyle \frac{C^2}{R_{}^4}}=0,`$ (9)
$`2m+{\displaystyle \frac{2KC_{}}{3}}+{\displaystyle \frac{2K^2R_{}^3}{9}}{\displaystyle \frac{4C_{}^2}{R_{}^3}}=0.`$ (10)
One can easily show, if $`C_{}`$ and $`R_{}`$ are critical, then the sign of

$$\beta 2C_{}+\frac{2}{3}KR_{}^3$$
(11)
is the same as the sign of $`C_{}`$.

There exists critical values of $`C_{}`$ that are positive ($`C_{}^+`$) or negative ($`C_{}`$). For definiteness we shall consider only the case when $`C(t=0)>C_{}`$, therefore the only limiting case we consider is that with $`CC_{}^+`$. (That choice corresponds to a foliation formed by leaves connecting two null infinities which moves forward in time - see a discussion in Sec IV in ). For simplicity we will drop the <sup>+</sup> suffix and $`C_{}`$ will mean a positive critical parameter. From the dynamical equation (6) follows that $`C`$ can only increase.

Next, let us introduce the notation that

$`ϵC_{}C`$ (12)
$`R_0R_{}+\delta .`$ (13)
where both $`\delta `$ and $`ϵ`$ are positive and small.

The equation $`p(R_0)=0`$ yields a nonlinear algebraic equation whose truncation gives

$$\delta ^2A+ϵ\beta =0.$$
(14)
Here $`A2R_{}^2+K^2R_{}^4`$. Eq. (14) is in fact the Lyapunov - Schmidt reduced equation constructed according to the standard rules . Therefore in the vicinity of the critical point we have

$$\delta =\sqrt{\frac{\beta ϵ}{A}}.$$
(15)
The function $`p`$ can be expressed in a form

$$\frac{pr}{2}=\sqrt{1\frac{R_0}{r}}\left[\frac{\kappa \delta }{R_0}+\frac{K^2}{9}(rR_0+r^22R_0^2)\frac{C^2}{R_0^4}(\frac{R_0}{r}+\frac{R_0^2}{r^2}+\frac{R_0^3}{r^3}3)\right]^{1/2}.$$
(16)
The insertion of (13), (15) and (16) into the equation (7) and the change of the integration variable to $`y=\sqrt{1\frac{R_0}{r}}`$, yield after a simple but tedious algebra

$$I(R_0)\sqrt{R_{}}_0^1𝑑y\frac{F_1}{\sqrt{\kappa \delta +y^2F_2}(\kappa \delta +y^2F_3)^2}.$$
(17)
Here the functions $`F_1,F_2`$ and $`F_3`$ are given by

$`F_1(y)`$ $`=`$ $`{\displaystyle \frac{K^2R_{}^3}{3(1y^2)^4}}+{\displaystyle \frac{6C_{}^2}{R_{}^3}}(1y^2)^2,`$ (18)
$`F_2(y)`$ $`=`$ $`{\displaystyle \frac{K^2R_{}^3(32y^2)}{9(1y^2)^2}}+{\displaystyle \frac{C_{}^2}{R_{}^3}}\left(64y^2+y^4\right),`$ (19)
$`F_3(y)`$ $`=`$ $`(33y^2+y^4)\left({\displaystyle \frac{2K^2R_{}^3}{9(1y^2)^3}}+{\displaystyle \frac{4C_{}^2}{R_{}^3}}\right),`$ (20)
while $`\kappa `$ reads

$$\kappa =\frac{2K^2R_{}^2}{3}+\frac{12C_{}^2}{R_{}^4}.$$
(21)
## E Limiting behaviour of the foliation.

The asymtotic behaviour of $`C`$ will be dominated by the $`1/\delta ^2`$ part of $`I(R_0)`$. As will be shown later, $`C`$ tends exponentially to $`C_{}`$; the attenuation factor in the exponent depends only on the leading term of $`I(R_0)`$. It is useful to define $`z=y/\sqrt{\kappa \delta }`$. Then one obtains $`I(R_0)=\frac{\sqrt{R_{}}}{\kappa ^2\delta ^2}\times I_d`$, where

$$I_d_0^{1/\sqrt{\kappa \delta }}𝑑z\frac{F_1(\sqrt{\kappa \delta }z)}{\sqrt{1+\frac{z^2}{R_{}}F_2(\sqrt{\kappa \delta }z)}(1+\frac{z^2}{R_{}}F_3(\sqrt{\kappa \delta }z))^2}.$$
(22)
One can split the integral $`_0^{1/\sqrt{\kappa \delta }}`$ into two parts: $`_0^{1/\sqrt{\kappa \delta }}=_0^{1/\sqrt{10^4\kappa \delta }}+_{1/\sqrt{10^4\kappa \delta }}^{1/\sqrt{\kappa \delta }}`$. It is easy to check that the contribution coming from the second integral goes to zero as $`\delta `$ approaches zero. Therefore $`F_1R_{}\kappa /2`$, $`F_2R_{}\kappa /2`$ and $`F_3R_{}\kappa `$. Thus the first integral (and also the integral $`I_d`$) is well approximated by

$`I`$ $`=`$ $`{\displaystyle \frac{\kappa R_{}}{2}}{\displaystyle _0^{\mathrm{}}}𝑑z{\displaystyle \frac{1}{\sqrt{1+\frac{R_{}\kappa }{2}z^2}(1+R_{}\kappa z^2)^2}}=`$ (24)
$`{\displaystyle \frac{\sqrt{\kappa R_{}}}{2}}{\displaystyle _0^{\mathrm{}}}𝑑z{\displaystyle \frac{1}{\sqrt{1+\frac{z^2}{2}}(1+z^2)^2}}.`$
The integral $`I_z=_0^{\mathrm{}}𝑑z\frac{1}{\sqrt{1+\frac{z^2}{2}}(1+z^2)^2}`$ can be explicitly evaluated and gives $`\sqrt{2}/2`$.

In summary, near the critical point we have

$`I(R_0)`$ $`=`$ $`{\displaystyle \frac{\sqrt{R_{}}}{\kappa ^2\delta ^2}}I`$ (26)
$`{\displaystyle \frac{\sqrt{2}}{4ϵ}}{\displaystyle \frac{R_{}A}{\kappa ^{3/2}|\beta |}}`$
The insertion of (13) and (26) into (6) yields

$$_tϵ=\mathrm{\Gamma }ϵ,$$
(27)
where

$$\mathrm{\Gamma }=\frac{|\beta |\kappa ^{3/2}}{2\sqrt{2}AR_{}}.$$
(28)
Eq. (27) immediately implies that $`ϵ`$ approaches 0 exponentially as

$$ϵ(t)=ϵ_0e^{t\mathrm{\Gamma }},$$
(29)
where $`ϵ_0`$ is an initial value of the parameter. Taking into account relations (13) and (15), one can conclude that the parameter $`C`$ and the minimal radius $`R_0`$ tend exponentially to their critical values $`C_{}`$ and $`R_{}`$, respectively.

Finally, collecting the above information and putting it into Eq. (8), we obtain the asymptotic behaviour of the lapse function near the critical point

$$N=N_0e^{t\mathrm{\Gamma }/2}.$$
(30)
This is exactly the same result as that obtained in ; in order to show equivalence, use expressions for the extrinsic curvatures (valid in the upper quadrant of the extended Schwarzschild geometry), which imply $`K=(2mR_{}3m)/\sqrt{2mR_{}^3R_{}^4}`$ and $`C_{}=(3mR_{}^3R_{}^4)/\sqrt{2mR_{}^3R_{}^4}`$. In the case of maximal slicing ($`K=0`$) the decay constant $`\mathrm{\Gamma }/2`$ is equal to $`4/(3\sqrt{6})`$, in agreement with the analytic derivation of  and close to the numerical result of  . The asymptotic behaviour of $`\gamma `$ and $`p`$ in a region close to the line $`R_0`$ is given by

$`\gamma `$ $`=`$ $`\gamma _0{\displaystyle \frac{e^{\frac{t\mathrm{\Gamma }}{4}}}{\sqrt{1\frac{R_0}{r}}}}`$ (31)
$`{\displaystyle \frac{pr}{2}}`$ $`=`$ $`p_0\sqrt{1{\displaystyle \frac{R_0}{r}}}e^{\frac{t\mathrm{\Gamma }}{4}}`$ (32)
Here $`\gamma _0`$ and $`p_0`$ are initial values of $`\gamma \sqrt{1\frac{R_{}}{r}}`$ and $`pr/(2\sqrt{1\frac{R_{}}{r}})`$, respectively. The four constants ($`ϵ_0,N_0,\gamma _0`$ and $`p_0`$) can be expressed in terms of one free parameter (say, $`ϵ_0`$) and $`A,\beta ,\kappa `$. Equations (3), (17), (29) - (32) and $`C=C_{}ϵ_0e^{t\mathrm{\Gamma }}`$ suffice to construct the metric (4) near the line $`R_0(t)`$ of minimal surfaces.

Acknowledgments. This work has been suported in part by the KBN grant 2 PO3B 006 23.