text:
# Dynamic exchange coupling and Gilbert damping in magnetic multilayers

## I Introduction

The spin transfer in ferromagnet–normal-metal hybrids causes a number of exciting phenomena. Equilibrium spin currents across a thin normal-metal (N) layer separating two ferromagnetic (F) films explain the oscillating RKKY-type exchange coupling Slonczewski (1993) as a function of the spacer thickness. This nonlocal exchange interaction may stabilize an antiparallel equilibrium configuration in magnetic multilayers Grünberg et al. (1986), which display a giant magnetoresistance Baibich et al. (1988). For thicker N spacers or tunneling barriers, this static ferromagnetic exchange vanishes, but the F layers can still be made to communicate by driving a dc current through the system. For example, depending on the direction of the current flow perpendicular to a layered F$``$N$``$F spin-valve structure, a torque can be exerted on the magnetizations. When one magnetization is fixed, this torque favors either parallel or antiparallel configurations, and may lead to a switching of the other magnetization Slonczewski (1996, 1999). Transport in the dynamic regime of moving magnetization directions has attracted relatively little attention, however, and is the subject of this paper.

We have recently launched the idea of spin pumping in weakly excited magnetic nanostructures. A moving ferromagnetic magnetization emits spins into adjacent conductors, exerting a relaxation torque and transferring an angular momentum out of the ferromagnet. The ejected spin angular momentum can either scatter back, relax in a nonmagnetic spacer, or be absorbed by a second ferromagnet. In the first case, the magnetization dynamics are not affected; in the second case, the magnetization motion is (nonlocally) damped by spin-flip scattering processes in the normal metal; and in the last case, the two ferromagnets become coupled by an exchange of itinerant spins, which may result in collective excitation modes. The spin-pumping concept has proven to be fruitful in the understanding of both the nonlocal damping mechanism Tserkovnyak et al. (2002a) and dynamic exchange interaction Heinrich et al. (2002) in hybrid F$``$N systems. The functionality of magnetic devices, such as magnetic random-access memories, is strongly affected by relaxation characteristics of the magnetic media. The nonlocal relaxation and coupling can be large and even dominant compared with other mechanisms in ultrathin films, and can be engineered by composition, geometry, and magnetic configuration of hybrid F$``$N systems. We hope that our theory will be useful as a designer tool, helping to realize future nanoscale magnetic devices with improved speed, operation threshold bias, and power consumption.

## II Enhanced damping in single F films

Consider an F$``$N bilayer as in Fig. 1. Without a voltage bias, no spin or charge currents flow when the magnetization of the ferromagnet is constant in time. When the magnetization direction starts precessing (as, e.g., under the influence of an applied magnetic field), a spin current $`𝐈_s^{\text{pump}}`$ is pumped out of the ferromagnet into the N layer Tserkovnyak et al. (2002a). When the ferromagnetic film is thicker than its transverse spin-coherence length $`\lambda _{\text{sc}}`$, $`d>\lambda _{\text{sc}}=\pi /|k_\text{F}^{}k_\text{F}^{}|`$, $`k_\text{F}^{()}`$ being the spin-dependent Fermi wave vectors, this current depends on the complex-valued parameter $`g^{}=g_r^{}+ig_i^{}`$, the interfacial “mixing” conductance Brataas et al. (2000); Waintal et al. (2000), by

$$𝐈_s^{\text{pump}}=\frac{\mathrm{}}{4\pi }\left(g_r^{}𝐦\times \frac{d𝐦}{dt}g_i^{}\frac{d𝐦}{dt}\right).$$
(1)
Here the time-dependent order parameter of the ferromagnet is a unit vector $`𝐦(t)`$, assuming a monodomain magnet with a spatially uniform magnetization at all times. A detailed derivation of Eq. (1) based on the scattering-matrix theory is given in Refs. Tserkovnyak et al. (2002a, b). Alternatively, this result can be derived in the framework of magnetoelectronic circuit theory Brataas et al. (2000, 2001) using only energy and angular-momentum conservation Tserkovnyak et al. (2002b). The dc conductance matrix is defined by

$$g^{\sigma \sigma ^{}}=\underset{mn}{}\left[\delta _{mn}r_{mn}^\sigma (r_{mn}^\sigma ^{})^{}\right].$$
(2)
Here $`r_{mn}^{}`$ ($`r_{mn}^{}`$) is the reflection coefficient off the F$``$N interface for spin-up (spin-down) electrons on the normal-metal side, which can be obtained by ab initio band-structure calculations Xia et al. (2002). $`g^{\sigma \sigma ^{}}`$ defined in Eq. (2) must be renormalized for highly transparent interfaces, as explained in Ref. Bauer et al. (2002), i.e., in the case of (globally) diffuse transport, the spurious Sharvin contributions have to be stripped off the real-valued conventional (Landauer-Büttiker) resistances for a definite spin $`\sigma `$:

$$\frac{1}{\stackrel{~}{g}^{\sigma \sigma }}=\frac{1}{g^{\sigma \sigma }}\frac{1}{2}\left(\frac{1}{N_N}+\frac{1}{N_{F\sigma }}\right),$$
(3)
while the inverse mixing conductance has to be corrected for the normal-metal Sharvin contribution only:

$$\frac{1}{\stackrel{~}{g}^{}}=\frac{1}{g^{}}\frac{1}{2N_N}.$$
(4)
The tildes in Eqs. (3) and (4) denote the renormalized conductances, which reduce to the bare values (2) when the numbers of transverse channels in the normal metal, $`N_N`$, and the ferromagnet, $`N_{F\sigma }`$, sufficiently exceed the contact conductances $`g^{\sigma \sigma ^{}}`$.

The total spin current, $`𝐈_s`$, across the F$``$N interface also has a backflow contribution, $`𝐈_s^{\text{back}}`$, in addition to the pumped current, $`𝐈_s^{\text{pump}}`$, see Fig. 1. The total spin transfer in the steady state,

$$𝐈_s=𝐈_s^{\text{pump}}𝐈_s^{\text{back}},$$
(5)
is determined self-consistently by the spin accumulation close to the interface. A finite value of the spin current, $`𝐈_s`$, would indicate the presence of a spin-sink mechanism in the normal metal by, e.g., spin-orbit coupling, when the spin angular momentum is transferred from the electron system to the lattice. In the case of small spin relaxation, the dynamically created spin accumulation in the normal metal may serve as a spin-battery device Brataas et al. (2002). In the opposite spin-flip scattering regime, the angular-momentum loss of the ferromagnet by $`𝐈_s`$ results in a Gilbert-type damping of the magnetic precession Tserkovnyak et al. (2002a).

The spin current out of the ferromagnet carries angular momentum perpendicular to the magnetization direction, corresponding to a torque $`𝝉=𝐈_s`$ on the ferromagnetic condensate Slonczewski (1996). Disregarding interfacial spin flips, this torque is entirely transferred to the magnetization, which is described by a generalized Landau-Lifshitz-Gilbert (LLG) equation Gilbert (1955); Slonczewski (1996)

$$\frac{d𝐦}{dt}=\gamma 𝐦\times 𝐇_{\text{eff}}+\alpha _0𝐦\times \frac{d𝐦}{dt}+\frac{\gamma }{M_sV}𝐈_s,$$
(6)
where $`\gamma `$ is the absolute value of the gyromagnetic ratio, $`\alpha _0`$ is the dimensionless intrinsic Gilbert-damping constant, $`M_s`$ is the saturation magnetization of the ferromagnet, and $`V`$ is its volume. Referring to Eq. (1) we see that the real part of the mixing conductance contributes to the damping just like the intrinsic bulk constant $`\alpha _0`$ which is thus smaller than the total Gilbert damping $`\alpha =\alpha _0+\alpha ^{}`$, whereas the imaginary part of the mixing conductance contributes like an effective field. The additional damping $`\alpha ^{}`$ is observable in, for example, ferromagnetic-resonance (FMR) spectra. Although not necessarily so for ferromagnetic insulators Huertas-Hernando et al. (2002), the mixing conductance for intermetallic F$``$N interfaces is to a good approximation real Xia et al. (2002) and therefore

$$\alpha ^{}=\kappa \frac{\gamma \mathrm{}\stackrel{~}{g}^{}}{4\pi M_sV}.$$
(7)
$`\kappa =1`$ corresponds to the perfect spin-sink model, when all pumped spins relax in the N layer and the backflow $`𝐈_s^{\text{back}}`$ vanishes. $`\kappa <1`$ corresponds to a finite backflow Tserkovnyak et al. (2002b).

The spin-sink capacity of clean normal-metal layers in contact with a ferromagnetic film, as in Fig. 1, is usually governed by spin-orbit scattering processes at impurities or defects. As the spin-flip probability $`ϵ=\tau _{\text{el}}/\tau _{\text{so}}`$ (defined in terms of the elastic scattering time $`\tau _{\text{el}}`$ and spin-orbit relaxation time $`\tau _{\text{so}}`$) rapidly increases with the atomic number $`Z`$, $`ϵZ^4`$ Abrikosov and Gor’kov (1962); Meservey and Tedrow (1978), we expect a larger spin-sink effect \[and therefore $`\kappa `$ in Eq. (7)\] for heavier metals (impurities as well as hosts). The extent of the hybridization of the conduction bands with $`p`$ or $`d`$ orbitals also plays an important role. In particular, clean noble metals, Cu, Ag, and Au, with predominantly $`s`$ conduction electrons are poor spin sinks with correspondingly small $`ϵ`$, but Pd and Pt, whose conduction electrons have significant $`d`$ character, have a high spin-orbit scattering rates and are efficient spin sinks. We note that heavy or magnetic impurities can turn an otherwise poor spin sink into a good one. The hierarchy of the Gilbert-damping enhancement has been measured for normal-metal buffers of Cu, Ta, Pd, and Pt (in order of increasing damping) Mizukami et al. (2001), in agreement with the aforementioned arguments Tserkovnyak et al. (2002b). In particular, the effect of spin pumping on the magnetization dynamics was shown to be negligible in the case of Cu, while Pt was a nearly perfect spin sink, resulting in a large increase of the Gilbert constant; this was also experimentally confirmed by Ingvarsson et al. Ingvarsson et al. (2002).

To further investigate the spin-sink effect of the normal-metal buffer, we next consider a more complex system, consisting of a double N layer attached to the ferromagnet, see Fig. 2. An interesting situation arises when the layer N1 is a bad spin sink, such as Cu, and N2 is a perfect spin sink, such as Pt Tserkovnyak et al. (2002b). Mizukami et al. Mizukami et al. (2002, 2002) experimentally studied the FMR line width in permalloy (Py)$``$Cu$``$Pt composites as a function of Cu (N1) width $`L`$. Next to $`L`$, there are three relevant length scales in the problem: the Fermi wave length, $`\lambda _F`$, the elastic scattering mean free path, $`\lambda _{\text{el}}`$, and the spin-diffusion length, $`\lambda _{\text{sd}}`$. In the coherent regime, $`L<\lambda _F`$, we do not expect large effects of a dusting N1 layer on the mixing conductance, which means that the damping enhancement (7) is approximately governed by $`\stackrel{~}{g}_{FN2}^{}`$ (with a quantum modulation when $`L\lambda _F`$). If $`L>\lambda _F`$, the spin transport across the N1 spacer may be described by the diffusion equation Tserkovnyak et al. (2002b), provided that either the spacer or the interface is disordered. In this regime, the spin backflow can be partitioned between the ferromagnet and layer N2, see Fig. 2. The relevant effective mixing conductance $`\stackrel{~}{g}_{\text{eff}}^{}`$ then has to account for scattering at both interfaces, F$``$N1 and N1$``$N2, as well as in the N1 spacer:

$$\frac{1}{\stackrel{~}{g}_{\text{eff}}^{}}=\frac{1}{\stackrel{~}{g}_{FN1}^{}}+R_{N1}+\frac{1}{\stackrel{~}{g}_{N1N2}},$$
(8)
where $`R_{N1}`$ is the resistance (per spin, in units of $`h/e^2`$) of the N1 layer, $`\stackrel{~}{g}_{N1N2}`$ is the one-spin conductance of the N1$``$N2 interface, and we assumed that $`L<\lambda _{\text{sd}}`$. If in addition $`L<\lambda _{\text{el}}`$, the bulk scattering, $`R_{N1}`$, can be disregarded and the total resistance is given simply by the sum of the interfacial contributions. $`\lambda _F`$ thus sets the length scale for a sharp drop in the Gilbert-damping enhancement, followed by an algebraic decay with the effective mixing conductance (8) for $`L>\lambda _F`$. Such a damping drop was reported in Ref. Mizukami et al. (2002) for Py$``$Cu$`(L)`$$``$Pt hybrids.

The regime $`L>\lambda _F`$ was studied in Refs. Mizukami et al. (2002, 2002), where a smooth algebraic decay of the damping enhancement was measured for $`L<\lambda _{\text{sd}}`$ followed by an exponential suppression for thicker Cu spacers, in excellent agreement with our theory Tserkovnyak et al. (2002b). Mizukami et al. Mizukami et al. (2002) also offered an explanation for their measurements, using a phenomenological theory due to Silsbee et al. Silsbee et al. (1979). The theoretically calculated damping profiles, $`\alpha (L)`$, as a function of $`L`$ in Refs. Tserkovnyak et al. (2002b) and Mizukami et al. (2002) are barely distinguishable. In our opinion, the phenomenological approach employed in the latter has severe limitations, however. In particular, the separation between localized and conduction electron spins Silsbee et al. (1979); Mizukami et al. (2002) is not justified for itinerant ferromagnets like the transition metals. Furthermore, this separation leads to wrong results for insulating ferromagnets which, we believe, can generate a spin current into adjacent nonmagnetic conductors in the same way as a conducting ferromagnet Tserkovnyak et al. (2002b), in contrast to the prediction of Ref. Silsbee et al. (1979). The good fit of the experimental $`\alpha (L)`$ in Ref. Mizukami et al. (2002) reflects the diffuse nature of transport in the copper layer rather than the spin-injection mechanism at the F$``$N interface. The magnitude of the spin current as calculated in Ref. Mizukami et al. (2002) is governed by a fitting parameter ($`\mathrm{\Gamma }`$ in their notation) which cannot be verified or predicted by ab initio calculations, unlike our spin-pumping expression (1).

The generality of our approach is demonstrated in the following by showing its validity for a qualitatively different system of two magnetic layers separated by a normal-metal film with negligible spin flip. We find that the same mechanism (namely the spin pumping) is responsible for the damping enhancement in coupled ferromagnets as well as the F$``$N structures discussed so far.

## III Dynamic coupling in spin valves

In F$``$N$``$F spin valves, the spin pumping causes qualitatively different effects in addition to the ones just described. Consider a system shown in Fig. 3. In the following we take $`L>\lambda _F`$, so that quantum coherence (in particular the static exchange coupling) can be disregarded, but $`L<\lambda _{\text{el}}`$, so that the spins of electrons are transferred between the two ferromagnets ballistically.

The total spin current pumped into the normal metal,

$$𝐈_s^{\text{pump}}=\frac{\mathrm{}}{4\pi }\left(\stackrel{~}{g}_1^{}𝐦_1\times \frac{d𝐦_1}{dt}+\stackrel{~}{g}_2^{}𝐦_2\times \frac{d𝐦_2}{dt}\right),$$
(9)
has contributions from both Fi$``$N interfaces, see Eq. (1). Here the small imaginary part of the mixing conductance Xia et al. (2002) is again disregarded. In a collinear configuration, the spin accumulation $`𝝁_s`$ induced by the spin pumping is perpendicular to the magnetizations. This simplifies the expression for the bias-driven spin current to Brataas et al. (2001); Tserkovnyak et al. (2002b); Brataas et al. (2002)

$$𝐈_s^{\text{back}}=\frac{\stackrel{~}{g}_1^{}+\stackrel{~}{g}_2^{}}{4\pi }𝝁_s.$$
(10)
A steady state is established when the two spin currents, Eqs. (9) and (10), cancel each other: $`𝐈_s^{\text{back}}=𝐈_s^{\text{pump}}`$. The spin accumulation is then given by $`𝝁_s=4\pi (\stackrel{~}{g}_1^{}+\stackrel{~}{g}_2^{})^1𝐈_s^{\text{pump}}`$. In general, $`𝐈_{s2}=𝐈_{s1}`$, by the conservation of angular momentum, but the net spin current $`𝐈_{sj}=𝐈_{sj}^{\text{pump}}𝐈_{sj}^{\text{back}}`$ across a single interface does not vanish. For the left ($`j=1`$) interface we find

$$𝐈_{s1}=\frac{\mathrm{}}{4\pi }\frac{\stackrel{~}{g}_1^{}\stackrel{~}{g}_2^{}}{\stackrel{~}{g}_1^{}+\stackrel{~}{g}_2^{}}\left(𝐦_1\times \frac{d𝐦_1}{dt}𝐦_2\times \frac{d𝐦_2}{dt}\right).$$
(11)
The normal-metal layer in our model scrambles the incoming spin current and divides it back over both ferromagnets. \[For very transparent interfaces the scrambling is only partial, but the present treatment is still adequate if we properly renormalize the conductance parameters as in Eqs. (3) and (4).\] The second ferromagnet can thus cause damping of the precession in the first magnetic film, and vice versa. A ferromagnetic layer can therefore serve as an efficient spin sink, just like a normal metal with strong spin-flip scattering.

We first discuss the implications of the spin transfer (11) when the second ferromagnet, F2, is stationary, $`\dot{𝐦}_2=0`$. In the FMR measurements, this is the case when one ferromagnet is in resonance, whereas the resonance frequency of the second ferromagnet is sufficiently different because of different magnetic anisotropies. The spin torque $`𝝉=𝐈_{s1}`$ on the first ferromagnet (F1) then has the form of the Gilbert damping when added as a source term to the LLG equation (6). In this case, the dynamic coupling of the ferromagnetic layers simply leads to an enhancement of the Gilbert-damping parameter with respect to its intrinsic value, exactly like in the case of single F films. The damping enhancement, following from Eq. (11), of the F1 layer therefore reads

$$\alpha ^{}=\frac{\gamma \mathrm{}\stackrel{~}{g}_1^{}\stackrel{~}{g}_2^{}}{4\pi M_sV(\stackrel{~}{g}_1^{}+\stackrel{~}{g}_2^{})}$$
(12)
(with $`\gamma `$, $`M_s`$, and $`V`$ of F1). Eq. (12) satisfactorily explains the increased viscous damping observed in Fe$``$Au$``$Fe spin valves Urban et al. (2001). The physical nature of the damping in single F films, Eq. (7), versus spin valves, Eq. (12), is quite clear. The angular momentum is first driven out of the ferromagnet via spin pumping. In the former case (single ferromagnet), the subsequent spin-orbit processes in the normal metal relax injected spins and, as a result, slow down the coherent motion of the ferromagnet. In the latter case (spin valve), the spin angular momentum is transferred from ferromagnet F1 and into F2. This spin current is absorbed by asynchronously driving the magnetization dynamics of F2. Therefore, both the normal-metal and ferromagnetic spin sinks act as an external brake to slow down the precession of the resonantly excited magnet. The spin-sink efficiency of the normal metal is characterized by the spin-flip probability $`ϵ`$ and its thickness $`L`$ Tserkovnyak et al. (2002b), which is lumped into the parameter $`\kappa `$ in Eq. (7). It follows from Eq. (12) that the efficiency of the adjacent ferromagnet F2, on the other hand, only depends on the F2$``$N mixing conductance (as long as the second ferromagnet is thicker than its transverse spin-coherence length $`\lambda _{\text{sc}}`$).

In the following we consider dynamics of the coupled F1$``$N$``$F2 system when both magnetizations are allowed to precess, i.e., when the ferromagnet resonance conditions are close to each other. By augmenting the LLG equation (6) for $`𝐦_1`$ by the spin current source term (11) and with small variables $`𝐮_i=𝐦_i𝐡_i`$, where $`|𝐮_i|1`$ and $`𝐮_i𝐡_j`$, we obtain the linearized expression (assuming a circular precession)

$$\dot{𝐮}_1=\omega _1𝐡_1\times 𝐮_1+\alpha _0𝐡_1\times \dot{𝐮}_1+\alpha _1^{}\left(𝐡_1\times \dot{𝐮}_1𝐡_2\times \dot{𝐮}_2\right),$$
(13)
where $`\omega _i`$ is the resonance frequency of ferromagnet Fi, $`𝐡_i`$ is the unit vector in the direction of its effective magnetic field, and $`\alpha _i^{}=\gamma _i\mathrm{}\stackrel{~}{g}_1^{}\stackrel{~}{g}_2^{}/[4\pi M_{si}V_i(\stackrel{~}{g}_1^{}+\stackrel{~}{g}_2^{})]`$, subscript $`i`$ labeling corresponding quantities of Fi. The dynamics of the magnetization direction $`𝐦_2`$ is obtained by exchanging subscripts $`12`$ in Eq. (13).

For a spin valve in the parallel configuration, $`𝐡_i=𝐡`$, and identical resonance frequencies, $`\omega _i=\omega `$, the quantity $`𝐮=𝐮_1/\alpha _1^{}+𝐮_2/\alpha _2^{}`$ (which, up to a scaling factor, is the rf component of the total angular momentum) is affected only by the intrinsic bulk damping, $`\dot{𝐮}=\omega 𝐡\times 𝐮+\alpha _0𝐡\times \dot{𝐮}`$, while the difference $`\mathrm{\Delta }𝐮=𝐮_1𝐮_2`$ relaxes according to $`\mathrm{\Delta }\dot{𝐮}=\omega 𝐡\times \mathrm{\Delta }𝐮+\alpha 𝐡\times \mathrm{\Delta }\dot{𝐮}`$ with the enhanced (Gilbert) damping constant $`\alpha =\alpha _0+\alpha _1^{}+\alpha _2^{}`$. The dynamic coupling in the antiparallel configuration as well as in the parallel configuration when the resonance frequencies have a large mismatch, $`\mathrm{\Delta }\omega (\alpha _0+\alpha _1^{}+\alpha _2^{})\omega `$, favors relaxation of each magnetization towards the equilibrium configuration. In this case, the two modes corresponding to excitations of either ferromagnet decouple, each having an enhanced damping parameter $`\alpha _i=\alpha _0+\alpha _i^{}`$. This picture explains in detail the FMR profiles measured on the Fe$``$Au$``$Fe spin valves, not only in both symmetric and very asymmetric limits discussed herein but also in the intermediate regime of closely matched but different resonance frequencies in the parallel alignment Heinrich et al. (2002).

## IV Conclusions

We presented a simple yet (for diffuse systems) complete theory of spin-current emission into adjacent conductors, which is generated by the time-dependent ferromagnetic order parameters. In the case of a single magnetic film in contact with a normal-metal layer, spin-flip processes in the latter viscously slow down the coherent precession of the ferromagnet, in accordance with the Gilbert phenomenology. When a magnetic layer is connected to other ferromagnets via Ohmic contacts, the latter can act as a nonlocal brake by relaxing the transverse angular momentum of the incident spins. Our framework formulated in terms of spin pumps and spin sinks successfully explains ferromagnetic resonance experiments on various multilayer systems.

In addition, this approach predicts a dynamic exchange coupling between ferromagnetic films connected by normal metals. Precessing magnetizations feel each other through a metallic spacer by exchanging nonequilibrium spin currents. Components of the ferromagnetic magnetizations that precess incoherently (i.e., out of phase) experience an additional damping obeying the Gilbert phenomenology. In symmetric spin valves and periodic superlattices, the dynamic exchange coupling causes a weakly damped, synchronized collective motion of the magnetic order parameters.

###### Acknowledgements.

The authors are grateful to B. I. Halperin and B. Heinrich for stimulating discussions. This work was supported in part by the DARPA Award No. MDA 972-01-1-0024, FOM, the NEDO International Joint Research Grant Program “Nano-magnetoelectronics”, and the NSF Grant No. DMR 99-81283.

Figure 1. Y. Tserkovnyak et al., J. Appl. Phys.

Figure 2. Y. Tserkovnyak et al., J. Appl. Phys.

Figure 3. Y. Tserkovnyak et al., J. Appl. Phys.


format: text

text:
# References

The CLEO collaboration has recently presented data indicating an observation of the decay cascade $`\mathrm{{\rm Y}}(3S)\gamma \chi _{bJ}^{}\gamma \omega \mathrm{{\rm Y}}`$, involving the hadronic transitions $`\chi _{bJ}^{}\omega \mathrm{{\rm Y}}`$. The kinematical analysis of the data favors dominance of transitions trough the $`\chi _{b1}^{}`$ resonance with the branching ratio $`B(\chi _{b1}^{}\omega \mathrm{{\rm Y}})=2.3\pm 0.4\%`$. It is the purpose of this letter to point out that within the standard description of the hadronic transitions between the levels of the $`b\overline{b}`$ system it would be rather troublesome to accommodate a dominance of the $`J=1`$ $`\chi _{b1}^{}`$ resonance over that with $`J=2`$, $`\chi _{b2}^{}`$. Namely, the expected ratio of the absolute decay rates, $`\mathrm{\Gamma }(\chi _{b2}^{}\omega \mathrm{{\rm Y}})/\mathrm{\Gamma }(\chi _{b1}^{}\omega \mathrm{{\rm Y}})`$, is given by the ratio of the $`S`$ wave phase space factors, which is approximately 1.4. Thus allowing for the relative rate of the initial radiative transitions from the $`\mathrm{{\rm Y}}(3S)`$ resonance and also for the ratio of the total decay rates of $`\chi _{b2}^{}`$ and $`\chi _{b1}^{}`$, one would estimate the relative contribution of these resonances to the observed rate of the cascade decays approximately as

$$r_{2/1}\frac{\mathrm{\Gamma }(\mathrm{{\rm Y}}(3S)\gamma \chi _{b2}^{}\gamma \omega \mathrm{{\rm Y}})}{\mathrm{\Gamma }(\mathrm{{\rm Y}}(3S)\gamma \chi _{b1}^{}\gamma \omega \mathrm{{\rm Y}})}1.3\pm 0.3,$$
(1)
where the indicated uncertainty arises mainly from the present knowledge of the ratio of the total decay rates of the $`\chi _{bJ}^{}`$ resonances.

Furthermore, in the nonrelativistic limit of heavy $`b`$ quarks, where the spin of the heavy quark effectively decouples, the spin structure of the decay amplitudes $`\chi _{b1}^{}\omega \mathrm{{\rm Y}}`$ and $`\chi _{b2}^{}\omega \mathrm{{\rm Y}}`$ is fully fixed, which leads to a distinct prediction for the angular correlation, different for each of these decays, between the observed decay products of $`\omega `$ and $`\mathrm{{\rm Y}}`$. Thus an additional information on the ratio $`r_{2/1}`$ can be found in the already existing data by studying the angular correlation between the direction determined by the leptons ($`e`$ or $`\mu `$) emerging from $`\mathrm{{\rm Y}}\mathrm{}^+\mathrm{}^{}`$ and the plane determined by the momenta of the charged pions from the decay $`\omega \pi ^+\pi ^{}\pi ^0`$.

In order to analyze the amplitudes of the decays $`\chi _{bJ}^{}\omega \mathrm{{\rm Y}}`$ one can recall that the hadronic transitions of the discussed type are viewed as proceeding through emission of a soft gluonic field in the heavy quarkonium transition, which field materializes as light hadron(s), the $`\omega `$ resonance in the case considered here. For nonrelativistic heavy $`b`$ quarks the dominant interaction with a soft gluon field is that with the chromoelectric component $`\stackrel{}{E}^a`$, since the interaction with the cromomagnetic field is suppressed by an inverse power of the heavy mass, $`m_b^1`$. It is important for the present discussion that the quantum numbers of the $`\omega `$ resonance ($`J^{PC}=1^{}`$) make possible the $`\omega `$ production by exclusively the chromoelectric gluon field i.e. without the need for its chromomagnetic component. The amplitude of such production can be described by a generic symmetric amplitude (in the rest frame of $`\omega `$)

$$\omega (\stackrel{}{ϵ})|d_{abc}f(q_1,q_2,q_3)\left(\stackrel{}{E}^a(q_1)\stackrel{}{E}^b(q_2)\right)\stackrel{}{E}^c(q_3)\delta ^4(q_1+q_2+q_3p)𝑑q_1𝑑q_2𝑑q_3|0=A_\omega \stackrel{}{ϵ}^{},$$
(2)
where $`\stackrel{}{ϵ}`$ stands for the polarization amplitude of $`\omega `$, $`p=(m_\omega ,\stackrel{}{0})`$ is its four-momentum, $`d_{abc}`$ are the symmetric SU(3) constants, and the form factor $`f(q_1,q_2,q_3)`$ is totally symmetric in its arguments and is determined by (presently poorly known) details of the dynamics of the transition between the heavy quarkonium levels.

In the limit where the heavy quarkonium is considered as a compact object and only its chromoelectric dipole interaction with $`\stackrel{}{E}^a`$ is retained, the form factor $`f`$ would be reduced to a point-like in the coordinate space. However it is known that such approximation does not work well for hadronic transitions from the $`\mathrm{{\rm Y}}(3S)`$ resonance, in particular for the decay $`\mathrm{{\rm Y}}(3S)\pi \pi \mathrm{{\rm Y}}`$ (for a discussion see e.g. the review , the latest data are presented in the talk), which behavior might be attributed to a rather high excitation energy of the $`\mathrm{{\rm Y}}(3S)`$. Since the $`\chi _{bJ}^{}`$ resonances are only slightly below the $`\mathrm{{\rm Y}}(3S)`$, an application of the leading terms in the multipole expansion to hadronic transitions from them can be somewhat questionable as well. Fortunately, the assumption of the applicability of this limit is not required for the considerations presented here. What is really important is that it is sufficient to consider the interaction of the heavy quarks only with the chromoelectric component of the gluon field. In the leading nonrelativistic order this interaction depends only on the coordinate variables of the heavy quarks and does not depend on their spin. In the same order the wave functions of the heavy quarkonium states can be factorized into the spatial part, $`\psi `$, and the spinor part $`\zeta `$. For the $`\mathrm{{\rm Y}}`$ resonance ($`{}_{}{}^{3}S_{1}^{}`$ state) the spatial part is a scalar $`\psi _0`$ and the spinor part is a vector, $`\stackrel{}{\zeta }_0`$, while for the $`\chi _{bJ}^{}`$ resonances ($`{}_{}{}^{3}P_{J}^{}`$ states) the spatial part transforms as a vector $`\stackrel{}{\psi }_1`$ and the spinor part is a vector too: $`\stackrel{}{\zeta }_1`$. Since, as discussed, the spin of the $`b`$ quarks decouples from the emission of $`\omega `$, the $`S`$ wave111It can be noticed that the fact that the discussed transitions are due to the dominant chromoelectric interaction and can proceed in the $`S`$ wave somewhat makes up for the ‘Zweig rule’ suppression, in a semi-quantitative agreement with the noticeable observed rate. amplitude of the transitions from the $`\chi _{bJ}^{}`$ states can be written as

$$A(\chi _{bJ}^{}\omega \mathrm{{\rm Y}})=Aϵ_i^{}\psi _0^{}\zeta _{0j}^{}P^{(J)}\psi _{1i}\zeta _{1j},$$
(3)
where $`P^{(J)}`$ is the projector of the direct product $`\psi _{1i}\zeta _{1j}`$ on the states with definite total spin $`J`$. Clearly, the higher partial wave amplitudes for the emission of $`\omega `$ can be neglected, given that the transitions from the $`\chi _{b1}^{}`$ and $`\chi _{b2}^{}`$ are just above the threshold: $`\mathrm{\Delta }_1M(\chi _{b1}^{})M(\mathrm{{\rm Y}})m_\omega =12.3\pm 0.6MeV`$, $`\mathrm{\Delta }_2M(\chi _{b2}^{})M(\mathrm{{\rm Y}})m_\omega =25.6\pm 0.5MeV`$ (while the $`\chi _{b0}^{}`$ resonance is well below the kinematical threshold for the discussed transition).

The fact that the dependence on the spin of the $`b\overline{b}`$ pair in eq.(3) factorizes out, implies that the ratio of the transition rates from the $`J=2`$ and $`J=1`$ states is given by the ratio of the $`S`$ wave phase space factors:

$$R_{2/1}\frac{\mathrm{\Gamma }(\chi _{b2}^{}\omega \mathrm{{\rm Y}})}{\mathrm{\Gamma }(\chi _{b1}^{}\omega \mathrm{{\rm Y}})}=\sqrt{\frac{\mathrm{\Delta }_2}{\mathrm{\Delta }_1}}1.4.$$
(4)
The ratio $`r_{2/1}`$ of the rates for the cascade processes can then be found as

$$r_{2/1}=\frac{\mathrm{\Gamma }(\mathrm{{\rm Y}}(3S)\gamma \chi _{b2}^{})}{\mathrm{\Gamma }(\mathrm{{\rm Y}}(3S)\gamma \chi _{b1}^{})}\frac{\mathrm{\Gamma }_{tot}(\chi _{b1})}{\mathrm{\Gamma }_{tot}(\chi _{b2})}R_{2/1}.$$
(5)
The ratio of the radiative decay rates from the $`\mathrm{{\rm Y}}(3S)`$ can be determined from the standard dipole-emission formula $`\mathrm{\Gamma }(\mathrm{{\rm Y}}(3S)\gamma \chi _{bJ}^{})(2J+1)k_\gamma ^3`$, which results in $`\mathrm{\Gamma }(\mathrm{{\rm Y}}(3S)\gamma \chi _{b2}^{})/\mathrm{\Gamma }(\mathrm{{\rm Y}}(3S)\gamma \chi _{b1}^{})1.1`$. The experimental value for this ratio is consistent with the theoretical number, but has a considerably larger uncertainty. A theoretical estimate of the ratio of the total decay rates of the $`\chi _{b2}^{}`$ and $`\chi _{b1}^{}`$ is plagued by the well known infrared sensitivity of the annihilation widths. This ratio however can be evaluated as an inverse of the relative branching ratia $`B(\chi _{bJ}\gamma \mathrm{{\rm Y}})`$222The data on the ratio of $`B(\chi _{bJ}\gamma \mathrm{{\rm Y}}(2S))`$ could be used as well. However the current experimental uncertainty in these data is larger than for the radiative transitions to $`\mathrm{{\rm Y}}(1S)`$. after a minor correction due to a slight difference in the factors $`k_\gamma ^3`$: $`\mathrm{\Gamma }_{tot}(\chi _{b1})/\mathrm{\Gamma }_{tot}(\chi _{b2})0.8\pm 0.15`$. Thus collecting the factors in eq.(5), one arrives at the estimate of the cascade rate ratio $`r_{2/1}`$ given in eq.(1).

Besides a slightly different kinematics in the discussed cascade transitions through the $`J=2`$ and $`J=1`$ resonances, their contribution can be additionally resolved by measuring the angular correlation between the decay products of the $`\mathrm{{\rm Y}}`$ and $`\omega `$ resonances. Indeed, these resonances are observed through their decays $`\mathrm{{\rm Y}}\mathrm{}^+\mathrm{}^{}`$ and $`\omega \pi ^+\pi ^{}\pi ^0`$. In the rest frame of the $`\mathrm{{\rm Y}}`$ resonance the two leptons are emitted ‘back to back’ thus defining one direction, and the other direction is the perpendicular to the plane defined by the charged pion momenta $`p_+`$ and $`p_{}`$. Defining the angle $`\theta `$ as that between these two directions, one can find the distribution of the rates over $`\mathrm{cos}\theta `$ for the cascade transitions through the $`\chi _{b2}^{}`$ and through the $`\chi _{b1}^{}`$, neglecting the motion of the $`\omega `$ in the rest frame of $`\mathrm{{\rm Y}}`$. (The small corrections due to this motion are proportional to $`\mathrm{\Delta }_{1,2}/m_\omega `$ and are parametrically of the same order as that due to the contribution of higher partial wave(s) in the amplitudes of $`\chi _{bJ}^{}\omega \mathrm{{\rm Y}}`$.) Since the structure of the decay amplitudes for $`\mathrm{{\rm Y}}`$ and $`\omega `$ is uniquely determined: $`A(\mathrm{{\rm Y}}\mathrm{}^+\mathrm{}^{})\stackrel{}{\zeta }_0(\overline{\mathrm{}}\stackrel{}{\gamma }\mathrm{})`$ and $`A(\omega \pi ^+\pi ^{}\pi ^0)\stackrel{}{ϵ}(\stackrel{}{p}_+\times \stackrel{}{p}_{})`$, the corresponding angular distributions are readily found from eq.(3) as

$$\frac{d\mathrm{\Gamma }(\chi _{b2}^{})}{d\mathrm{cos}\theta }\left(1\frac{1}{7}\mathrm{cos}^2\theta \right),\frac{d\mathrm{\Gamma }(\chi _{b1}^{})}{d\mathrm{cos}\theta }\left(1+\mathrm{cos}^2\theta \right).$$
(6)
The correlation coefficient $`a`$ in the overall angular distribution in the cascade $`\mathrm{\Gamma }(\mathrm{{\rm Y}}(3S)\gamma \chi _{bJ}^{}\gamma \omega \mathrm{{\rm Y}})(1+a\mathrm{cos}^2\theta )`$ with the two branches, corresponding to $`J=2`$ and $`J=1`$, mixed in the proportion determined by $`r_{2/1}`$, is therefore found as

$$a=\frac{5r_{2/1}}{5+7r_{2/1}}.$$
(7)
The expected value of the ratio $`r_{2/1}`$ in eq.(1) corresponds to $`a0.26\pm 0.06`$, while in the limit of a strong dominance of the cascade through $`\chi _{b1}^{}`$, i.e. when $`r_{2/1}0`$, the correlation coefficient obviously equals one. Hopefully, the already available data are sufficient for resolving between these two values.

I thank Jon Urheim, who brought to my attention the data, for enlightening discussion of the experimental situation and I gratefully acknowledge a useful discussion of the theoretical aspects with Arkady Vainshtein. This work is supported in part by the DOE grant DE-FG02-94ER40823.


text:
# Three-dimensional quasi-Tonks gas in a harmonic trap

## Abstract

We analyze the macroscopic dynamics of a Bose gas in a harmonic trap with a superimposed two-dimensional optical lattice, assuming a weak coupling between different lattice sites. We consider the situation in which the local chemical potential at each lattice site can be considered as that provided by the Lieb-Liniger solution. Due to the weak coupling between sites and the form of the chemical potential, the three-dimensional ground-state density profile and the excitation spectrum acquire remarkable properties different from both 1D and 3D gases. We call this system a quasi-Tonks gas. We discuss the range of applicability of this regime, as well as realistic experimental situations where it can be observed.

Recent developments in low-dimensional trapping 1D , loading of optical lattices Latt ; Greiner , and external modification of the interparticle interactions Feshbach , open fascinating perspectives towards the achievement of strongly-interacting systems in cold atomic gases. In this sense, a remarkable example is provided by the recent observation of the superfluid to Mott-Insulator (MI) transition in cold bosonic gases in optical lattices Greiner . Other strongly-interacting systems in cold atomic gases have been recently considered theoretically, as the case of large scattering length biga , rapidly rotating Bose gases belen , and one-dimensional systems, both bosonic Olshanii2 ; Chiara ; Letter ; Gangardt ; Giorgini  and fermionic Recati .

One-dimensional systems differ significantly when compared to 3D ones. Most remarkably, in 1D gases the interactions become relatively more important the more dilute the gas is. In particular, for very low densities or a large value of the $`s`$-wave scattering length $`a`$, the system enters the so-called Tonks-Girardeau (TG) regime, in which the bosons acquire fermionic properties Girardeau60 ; Wright . For finite interactions, the problem of delta-interacting 1D homogeneous Bose gases was exactly solved by Lieb and Liniger (LL) by means of Bethe Ansatz LiebLiniger . A combination of the LL solution, the local density approximation, and hydrodynamic equations have been recently employed in the analysis of different properties of trapped 1D Bose gases, as the ground-state density Olshanii2 , the excitation spectrum Chiara , the non self-similar expansion in a 1D guide Letter , and the correlation and coherence properties Gangardt . The latter has also been recently analyzed by means of quantum Monte Carlo methods Giorgini .

The accomplishment of strongly-interacting 1D gases requires tight transversal trapping, low atom numbers, and possibly the increase of $`a`$ via Feshbach resonances Olshanii2 ; Olshanii1 ; Petrov1D . In this sense, 2D optical lattices are specially favorable, since the on-site transversal confinement can be very strong. For sufficiently intense lattices, the tunneling rate could be made small enough to ensure that each lattice site behaves as an independent 1D experiment. This paper analyzes the situation in which the tunneling is non negligible, but sufficiently weak to assume the 1D character of the local chemical potential. In this regime, which we call the quasi-Tonks regime, we predict that the combination of tunneling and 1D local chemical potential results in important modifications of the 3D ground-state density and the excitations.

In the following we consider a gas of $`N`$ bosons in a cylindrically-symmetric harmonic trap, with axial (radial) frequency $`\omega _z`$ ($`\omega _{}`$). We assume a superimposed periodic potential provided by an optical lattice of the form $`V_\mathrm{l}(x)+V_\mathrm{l}(y)=V_0(\mathrm{sin}^2qx+\mathrm{sin}^2qy)`$, where $`q=2\pi /\lambda `$ is related to the laser wavelength $`\lambda `$, and $`V_0`$ is the lattice amplitude. The lattice has a periodicity $`d=\pi /q=\lambda /2`$. From a macroscopic point of view it is convenient to introduce the effective mass, $`m^{}`$, to describe the effects of the lattice on the dynamics. The value of $`m^{}`$ can be obtained from $`^2E_0/k^2=\mathrm{}^2/m^{}`$, where $`E_0(k)`$ is the dispersion law corresponding to the lowest energy band. In tight-binding regime, $`m^{}`$ can be related with the tunneling rate $`J`$, as $`m/m^{}=\pi ^2J/E_R`$, where $`E_R=\mathrm{}^2q^2/2m`$ is the recoil energy, $`J=w_i(x)(\mathrm{}^2_x^2/2m+V_\mathrm{l}(x))w_{i+1}(x)𝑑x`$, and $`\{w_i\}`$ are the Wannier functions for the lowest band.

We first consider the case of a single isolated lattice site, which can be well approximated as a dilute gas of $`N_t`$ bosons confined in a very elongated harmonic trap with radial and axial frequencies $`\omega _\rho `$ and $`\omega _z`$ ($`\omega _\rho \omega _z`$). The transversal confinement is strong enough to fulfill that the interaction energy per particle is smaller than the zero-point energy $`\mathrm{}\omega _\rho `$ of the transversal trap. In this way, the transversal dynamics is effectively “frozen” and the system can be considered 1D. The Hamiltonian describing the physics of the delta-interacting 1D gas is

$$\widehat{H}_{1\mathrm{D}}=\widehat{H}_{1\mathrm{D}}^0+\underset{j=i}{\overset{N_t}{}}\frac{m\omega _z^2z_i^2}{2}$$
(1)
where

$$\widehat{H}_{1\mathrm{D}}^0=\frac{\mathrm{}^2}{2m}\underset{j=1}{\overset{N_t}{}}\frac{^2}{z_j^2}+g_{1\mathrm{D}}\underset{i=1}{\overset{N_t1}{}}\underset{j=i+1}{\overset{N_t}{}}\delta \left(z_iz_j\right)$$
(2)
is the Hamiltonian in absence of the harmonic trap, and $`g_{1\mathrm{D}}=2\mathrm{}^2/ma_{1\mathrm{D}}`$. The scattering problem under one-dimensional constraints was analyzed in detail by Olshanii Olshanii1 , and it is characterized by the 1D scattering length $`a_{1\mathrm{D}}=(a_\rho ^2/2a)[1𝒞(a/a_\rho )]`$, with $`a`$ the 3D scattering length, $`a_\rho =\sqrt{2\mathrm{}/m\omega _\rho }`$ the oscillator length in the radial direction, and $`𝒞=1.4603\mathrm{}`$. For the thermodynamic limit, a 1D gas at zero temperature with linear density $`n`$, is characterized by an energy per particle, which can be obtained from the corresponding LL integral equations LiebLiniger . Assuming that the density variates sufficiently slowly, at each point $`z`$ the gas can be considered in local equilibrium, and the local energy per particle is provided by the LL equations for the density $`n(z)`$. We call this approach Local Lieb-Liniger (LLL) model.

Although, strictly speaking, the LLL approach is only valid for 1D systems, for sufficiently low $`J`$, the virtual processes involving transitions into nearest sites, should just result in a small correction to the 1D scattering of the order of $`J/\mu `$, where $`\mu `$ is the chemical potential. Therefore, if $`J/\mu 1`$, the actual local chemical potential can be well approximated by that obtained from the LLL approach. If the tunneling becomes very small, even for large values of $`N_t`$ the system could enter into the MI regime. However, this regime typically demands an extremely small value of $`J/\mu `$. Summarizing, in what we called the quasi-Tonks regime, $`1/N_tJ/\mu 1`$, the LLL chemical potential can be employed, and at the same time the tunneling cannot be neglected.

Assuming a sufficiently slow variation of the gas density, the dynamics can be well described by means of the corresponding macroscopic hydrodynamic equations meret :

$`{\displaystyle \frac{n}{t}}+(n𝐯)=0`$ (3)
$`m^{}{\displaystyle \frac{v_x}{t}}+{\displaystyle \frac{}{x}}(K+V+U)=0`$ (4)
$`m^{}{\displaystyle \frac{v_y}{t}}+{\displaystyle \frac{}{y}}(K+V+U)=0`$ (5)
$`m{\displaystyle \frac{v_z}{t}}+{\displaystyle \frac{}{z}}(K+V+U)=0`$ (6)
where

$`K`$ $`=`$ $`{\displaystyle \frac{m^{}}{2}}v_x^2+{\displaystyle \frac{m^{}}{2}}v_y^2+{\displaystyle \frac{m}{2}}v_z^2`$ (7)
$`V`$ $`=`$ $`{\displaystyle \frac{m}{2}}\omega _{}^2(x^2+y^2)+{\displaystyle \frac{m}{2}}\omega _z^2z^2`$ (8)
$`U`$ $`=`$ $`\mu (n)=\mu _{1D}(d^2n)`$ (9)
where $`n(x,y,z)`$ is the macroscopic 3D density and $`𝐯=(v_x,v_y,v_z)`$ is the velocity field. Note that all the information about the lattice is contained in $`m^{}`$, and in the local chemical potential $`\mu (n)`$. The use of the hydrodynamic equations is justified by the superfluidity of the quasi-1D tubes Sonin  and by the avoidance of the MI phase. In the equation (9) the local chemical potential of the 3D system is related to the 1D chemical potential, $`\mu _{1D}`$, provided by the LLL approach. In order to obtain $`\mu _{1D}`$ the value of $`a_{1D}`$ must be obtained for each site. The latter demands the knowledge of the effective oscillator length $`a_\rho `$ associated with the on-site radial confinement. Since the latter can be very large, a Gaussian Ansatz can be assumed for the on-site radial wave function. In this way, we can calculate $`a_\rho `$ by minimizing the radial energy

$`{\displaystyle \frac{E(a_\rho )}{E_R}}={\displaystyle \frac{1}{q^2a_\rho ^2}}+{\displaystyle \frac{V_0}{E_R}}\left(1e^{q^2a_\rho ^2}\right).`$ (10)
In Eq. (10) we have neglected the interaction energy, which for $`V_0/E_R>3`$ does not provide any significant contribution. Once known the value of $`a_{1D}`$, we evaluate $`\mu _{1D}(n)`$ by numerically solving the corresponding LL integral equations.

By imposing the equilibrium conditions $`n/t=v_i=0`$, we obtain the equation of state

$$\mu [n_0(x,y,z)]=\mu _TV(x,y,z),$$
(11)
where $`n_0`$ is the equilibrium density, and $`\mu _T`$ is the chemical potential of the system. Inverting (11), one obtains the expression for $`n_0`$. Imposing the normalization to the number of particles $`N`$ we obtain:

$`A={\displaystyle \frac{|a_{1D}|^4d^2N}{a_x^2a_{}^4}}=4\pi \xi ^3{\displaystyle _0^1}t^2\stackrel{~}{\mu }_{1D}^1[\xi ^2(1t^2)]𝑑t,`$ (12)
where $`\xi =2\mu m|a_{1D}|^2/\mathrm{}^2`$. From Eq. (12) we observe that, similarly to the strict 1D case Olshanii2 ; Chiara ; Gangardt , the ground state of the system is completely characterized by a single parameter $`A`$. For $`A1`$ the mean-field (MF) regime is retrieved, whereas the TG is found for $`A1`$. In Eq. (12), $`\stackrel{~}{\mu }_{1D}(x)=2m|a_{1D}|^2\mu _{1D}/\mathrm{}^2`$, is a function of the dimensionless density $`\eta =|a_{1D}|d^2n`$. Fig. 1 shows the density as a function of the dimensionless radius $`\stackrel{~}{r}`$, where $`\stackrel{~}{r}^2=_i(r_i/R_i)^2`$, with $`R_i`$ the Thomas-Fermi radii. In order to compare the different cases, we re-scale the density with respect to the central one. Nevertheless, we stress that the values of $`R_i`$ and the central density depend on the parameters, and in particular on the interaction regime.

From the previous discussion, it becomes clear that the 1D properties of the local chemical potential are transferred to the 3D stationary shape of the cloud. In particular, for small values of $`A`$, the transversal density profile is not an inverted parabola, as expected for the case of a BEC in a harmonic trap. We stress once more, that this remarkable effect is just possible in the presence of the 2D lattice, which guarantees, under the above discussed conditions, a LLL chemical potential.

In the following we calculate the frequencies for the lowest excitations of a quasi-Tonks gas. In particular, we evaluate the breathing and quadrupole modes, by using the Ansatz $`n=n_0\left(r_i/b_i\right)/_jb_j`$, $`v_i=(\dot{b}_i/b_i)r_i`$, which constitutes an exact solution of the continuity equation (3). Multiplying the Euler equations by $`nx_i`$ and integrating, one obtains

$`\ddot{b_x}+\stackrel{~}{\omega }_{}^2b_x+{\displaystyle \frac{\stackrel{~}{\omega }_{}^2}{b_x}}F({\displaystyle \underset{j=x,y,z}{}}b_j)=0`$ (13)
$`\ddot{b_y}+\stackrel{~}{\omega }_{}^2b_y+{\displaystyle \frac{\stackrel{~}{\omega }_{}^2}{b_y}}F({\displaystyle \underset{j=x,y,z}{}}b_j)=0`$ (14)
$`\ddot{b_z}+\omega _z^2b_z+{\displaystyle \frac{\omega _z^2}{b_z}}F({\displaystyle \underset{j=x,y,z}{}}b_j)=0,`$ (15)
where we have employed the effective frequencies $`\stackrel{~}{\omega }_{}^2=\omega _{}^2m/m^{}`$ and defined the function

$$F(\rho )=\frac{1}{m\omega _i^2Nx_i^2_0}n_0x_i\frac{}{x_i}\mu (\rho ^1n_0)d^3r.$$
(16)
Using the previously calculated equilibrium solution, one obtains that $`F(\rho )`$ is independent of the particular choice of the coordinate $`i`$, and that $`F(1)=1`$. Note that the function $`F`$ only depends on $`A`$, which, as previously discussed, completely characterizes the static properties of the system. Linearizing around the equilibrium solution $`b_i=1`$, one obtains the frequencies of the breathing mode ($`\omega _B`$), the $`m=0`$ quadrupole mode ($`\omega _{Q1}`$) and the $`m=2`$ quadrupole one ($`\omega _{Q2}`$):

$`\omega _B^2`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(\stackrel{~}{\omega }_{}^2(2+2f)+\omega _z^2(2+f)+\sqrt{[\stackrel{~}{\omega }_{}^2(2+2f)+\omega _z^2(2+f)]^28\stackrel{~}{\omega }_{}^2\omega _z^2(2+3f)}\right)`$ (17)
$`\omega _{Q1}^2`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(\stackrel{~}{\omega }_{}^2(2+2f)+\omega _z^2(2+f)\sqrt{[\stackrel{~}{\omega }_{}^2(2+2f)+\omega _z^2(2+f)]^28\stackrel{~}{\omega }_{}^2\omega _z^2(2+3f)}\right)`$ (18)
$`\omega _{Q2}^2`$ $`=`$ $`2\stackrel{~}{\omega }_{}^2`$ (19)
In the expression for the frequencies we employ the parameter $`f=F^{}(1)`$, which is an universal function of the parameter $`A`$. We show this dependence in Fig. 2. The parameter $`f`$ ranges from $`2`$, at the TG limit, to $`1`$, for the MF regime. In this two limiting cases the function $`F`$ can be obtained analytically, being $`F=\rho ^2`$ (TG) and $`F=\rho ^1`$ (MF).

The quasi Tonks regime is not only achievable for realistic conditions, but, actually, it is expected to be the case for typical parameters in ongoing experiments, as those of Fig. 3, which shows the dependence of the frequencies (17-19) on the lattice amplitude $`V_0/E_R`$. We stress that the system enters the MI regime for large values of $`V_0/E_R`$ (around $`40`$ for the case of Fig. 3). If this is the case, as discussed above, our macroscopic hydrodynamic approach fails, and the system becomes a set of independent 1D gases. On the other hand, for decreasing values of $`V_0/E_R`$ the system abandons the quasi-Tonks regime (which for the case of Fig. 3 occurs at $`V_0/E_R15`$), and the condition $`J/\mu 1`$ is reached. This constitutes an additional intermediate cross-dimensional regime, in which the gas is 3D, but the local chemical potential is not that expected for a 3D Bose gas. For this regime, contrary to the case of the quasi-Tonks gas, the local chemical potential cannot be approximated by the LLL chemical potential. As a consequence of that, the challenging problem of describing the $`J/\mu 1`$ regime demands a completely different approach. In this sense, especially interesting could be to employ the analogies to the problem of coupled Luttinger liquids Luttinger .

The existence of the quasi Tonks regime can be easily revealed in experiments, by either observing the size and the form of the stationary density profile, or monitoring the collective excitations. For increasing values of $`V_0/E_R`$, the frequency of the breathing mode approaches that of the lowest compressional mode in 1D gases Chiara , whereas the quadrupolar frequencies tend to zero. For intermediate values of the lattice potential, the lowest excitations are significantly different than in 1D, for the same value of the dimensionless density $`\eta `$. In particular, more than one excitation mode is expected. We stress that the modes (17-19) also differ significantly from the expected results for a BEC in a lattice meret . The latter are recovered from our formalism for large values of $`A`$ (MF regime). We must point, however, that the local correlation properties are preserved in the quasi Tonks regime. In this sense, if the parameter $`A`$ becomes sufficiently small the lifetime of the gas can be very significantly enlarged, due to the reduction of the two- and three-body losses Gangardt .

In this Letter, we have considered what we call the quasi-Tonks regime, in which a gas confined in a 2D optical lattice, can present significant tunneling, and at the same time maintain the local chemical potential obtained for each site using the LLL model. We have discussed that this regime should be the one encountered under current typical experimental conditions. In this cross-dimensional regime, a 3D cloud acquires 1D properties. In particular, we have shown that both the ground-state properties and the excitation spectrum are significantly different than the expected results for a harmonically confined BEC.

We thank M. Lewenstein, P. Öhberg, L. P. Pitaevskii, G. V. Shlyapnikov, S. Stringari, and the experimental Quantum Optics group of the University of Hannover for stimulating discussions. We acknowledge support from Deutsche Forschungsgemeinschaft SFB 407 and SPP1116, the RTN Cold Quantum Gases, IST Program EQUIP, ESF PESC BEC2000+, EPSRC, and the Ministero dell’Istruzione, dell’Università e della Ricerca (MIUR). The authors wish to thank the Alexander von Humboldt Foundation, the Federal Ministry of Education and Research and the ZIP Programme of the German Government.


text:
# Comment on “Estimating 1/𝑓^𝛼 scaling exponent from short time series”

## Abstract

A recently proposed method (O. Miramontes, P. Rohani, Physica D 166 (2002) 147) for estimating the scaling exponent in very short time series may give wrong results, especially in case of undersampled data.

Estimating the scaling exponent in the power spectrum of a signal is one of the most important tasks in the analysis of real–world temporal sequences. The point of time series analysis is not making a statement about the signal in question, but getting insight of the (unknown) mechanism responsible for generating the signal. One may say that the time series is usually not very interesting, but its generating mechanism is. The scaling exponent is one of the quantities that are calculated in order to get a better understanding of this mechanism. Many techniques have been developed for this job; see Ref.  for a recent review. It is particularly difficult to estimate the scaling exponent if the time series is very short, which is frequent for biological, medical and ecological data. Indeed, it is not clear whether the notion of the scaling exponent has any well–defined meaning in this case.

Short time series are usually characterized by a large sampling time, $`T`$. This is easy to understand: The series is short because it takes long time to gather the data, and it takes long time to gather the data because the sampling time is large and cannot be easily shortened. In a time series only the frequencies up to the Nyquist frequency $`f_\mathrm{N}=1/(2T)`$ are present . Any possible power scaling that is discovered within the series is necessarily restricted to the Nyquist interval $`[0,f_\mathrm{N}]`$. This is all the data can tell — there are no scientific grounds to claim that the scaling extends beyond the Nyquist frequency. Because $`T`$ is large, $`f_\mathrm{N}`$ is small. Therefore, attributing the power scaling, if one is detected in such a series, to noise is misleading at best, as in a true noise the power scaling extends to many (theoretically to infinitely many) orders of magnitude.

When dealing with a short time series, one is also faced with another, even more acute problem: the Nyquist interval is sparsely covered. There is little one can tell about the power density at frequencies that are not multiplies of $`1/(2NT)`$, where $`N`$ is the lenght of the series, and it is not clear whether the very notion of a power scaling can be used in such situations. Any attempts to interpolate between frequencies directly accessible form the data merely scale down the behaviour observed for these frequencies to those that are not directly accessible. This approach may work only if one knows a priori that the process used to generate the data is governed by a power law.

A method published recently in Ref.  suffers from all these deficiencies. The method, named “multiple segmenting method” (MSM), consists of two steps: First, one calculates scaling exponents for each segment of length $`2^s`$, $`s3`$, and calculates the average scaling exponents for each $`s`$; denote these averages by $`g(n)`$ with $`n=2^s`$. In this step pseudo–replicates of the original data are created. These pseudo-replicates have the same short–time, or large–frequency, correlations as the original series, and the process of replicating increases their statistical significance. However, any accidental correlations among the data are enhanced as well. Second, the relation

$$g(n)=a+\frac{b}{\sqrt{n}}$$
(1)
is fitted to the averages calculated in the first step, and the scaling exponent for the whole series is calculated from (1) by using results of the fit and $`n=N`$, where $`N`$ is the actual length of the time series. The equation (1) has not been derived in any rigorous way, only guessed from the example data.

The authors of Ref.  claim that MSM gives correct results for time series as short as $`N=47`$ terms, but offer little proof to this claim. They present a number of examples with artificially created time series of moderate length $`N=400`$, almost an order of magnitude larger than 47. The time series used had, by construction, built–in power scaling, and so did all their subseries. In these cases the MSM method discovered the scaling exponents that we beforehand knew were present in the data. Two out of four real–world examples presented had time series even longer ($`N=2000`$ and $`N=1024`$, respectively), where the power–type behaviour could be detected and the corresponding scaling exponents determined by conventional methods. There were only two really short ($`N=47`$) time series of unknown properties discussed. The MSM method, when applied to these two series, gave some numbers as the scaling exponents, but it is not clear whether these numbers are at all meaningful.

In what follows we present examples of short time series where the MSM method, when blindfoldedly applied, gives clearly misleading results.

As a first example, we take a quasiperiodic function

$$x(t)=\underset{i=1}{\overset{3}{}}a_i\mathrm{sin}(\omega _it),$$
(2)
with $`a_1=0.8`$, $`\omega _1=5`$, $`a_2=0.6`$, $`\omega _2=20\sqrt{3}`$, $`a_3=0.8`$, $`\omega _3=30\sqrt{3}`$, sampled with a time step $`T=7.9/460.17174`$ to give $`N=47`$. This time series is clearly undersampled, and many of the high–frequency features of the function (2) are lost, cf. Fig 1. When we apply the MSM method to this time series, we obtain a scaling exponent of about $`0.98`$, cf. Fig. 2 and Table 1. This result, if taken at its face value, indicates that the time series is governed by a flicker (or pink) noise, which is obviously untrue: the microscopic mechanism responsible for this time series is a regular, quasiperiodic function, albeit undersampled.

We used a realization of the Markovian AR(1) process

$$X_{n+1}=\beta X_n+(1\beta )\eta _n,$$
(3)
with $`\beta =0.69`$ like in , to generate two next example time series. $`\eta _n`$ is a Gaussian white noise, generated numerically by means of an algorithm published in . We generated 256 terms of the series (3), and then constructed our second example time series by taking every fifth term of the original series up to $`N=47`$ ($`X_5`$, $`X_{10}`$,…,$`X_{235}`$). Finally, the first consecutive 47 terms of the original series ($`X_1`$, $`X_2`$,…,$`X_{47}`$) were used as the third example. The output of the MSM method is presented in Figs. 4 and 5 and in Table 1. Note that the function (1) does not seem to be a reasonable fit to these data, and if this fit fails, MSM does not offer any other means to systematically determine the scaling exponent for the whole series. For the second example the method indicates the scaling exponent in the range $`[0.33,0.16]`$, i.e. it recognizes the time series as being close to the white noise, which is again not true. On the other hand, for the third example the method predicts the scaling exponent in the range $`[1.52,1.05]`$, which is close (especially the lowest value) to a rough estimate for the scaling exponent for the whole 256 term time series. Note that this time series was, by construction, characterized by a power–type behaviour. The relative success of MSM in this particular case strengthens our point that if a power law is present in the data, scaling it down to frequencies not observed in the data may not hurt much, but a short time series does not tell whether such a power law is present.

One may argue that the MSM method failed for the undersampled time series only. However, undersampling can be quite common in practice. When dealing with real–world time series generated by a mechanism whose details remain unknown, we never know whether we choose a correct sampling time. For instance, in the bacterial population data analysed in Ref.  the sampling time was as long as a week, while the true dynamics of the population might have been governed by processes with a much lower characteristic time, like daily changes in sunlight or temperature. Similarly, in many medical studies levels of drugs in subjects are determined on a daily basis, while physiological processes responsible for the drug absorption and decay usually proceed much faster. There are many other situations in which data are collected not in a controlled laboratory environment, but come out of real–world processes, where adjusting the sampling time may be impossible for various “technical” reasons. In case of undersamplig, the power spectrum beyond the Nyquist frequency is aliased into the Nyquist interval, distorting the available power spectrum . Examples presented above show that the MSM method, when applied to such data, may give wildly wrong results.

One needs to know some important details of the underlying process in order to estimate its parameters from a very short time series. Otherwise idiosyncrasies of the data are mistakenly taken for true correlations and are enhanced in the process of pseudo–replicating the time series. This may lead to a failure of the MSM method. This method appears to work tolerably well for long time series, but in this case it is inferior to conventional techniques due to its large computational costs. In our opinion, the usefulness of the MSM method is thus limited to the rare case when one knows for sure that the process used to generate the time series is governed by a power law with an unknown exponent, but the time series available is very short and cannot be easily extended and only a rough estimate of the scaling exponent is required. In other situations the MSM method is not recommended.

I would like to thank Dr. Adam Kleczkowski from Cambridge University for his helpful comments.