text:
# Recent Progress in Mathematical Analysis of Vortex Sheets

## 1. Introduction

Vortex dynamics is of fundamental importance for a wide variety of concrete physical problems, such as lift of airfoils, mixing of fluids, separation of boundary layers, and generation of sounds. In mathematical analysis, one often neglects surface tension and viscosity, when they are small in the real physical problem. This necessitates justifying such simplifications.

In this paper, we consider the motion of the interface separating two domains of the same fluid in $`R^2`$ that moves with different velocity along the tangential direction of the interface. We assume that the fluids occupying the two domains separated by the interface are of constant densities that are equal, are inviscid, incompressible and irrotational. We also assume that the surface tension is zero, and there is no external forces. The interface in the aforementioned fluid motion is a so called vortex sheet. We want to study the following problem:

Given a vortex sheet initial data, is there a unique solution to this problem?

In general, there are two approaches to the aforementioned problem. One is to solve the initial value problem of the incompressible Euler equation in $`R^2`$:

$$\{\begin{array}{cc}v_t+vv+p=0\text{div}v=0,\hfill & \\ v(x,y,0)=v_0(x,y)\hfill & \end{array}(x,y)R^2,t0$$
$`1`$
where the initial incompressible velocity $`v_0L_{loc}^2(R^2)`$, in which the vorticity $`\omega _0=\text{curl}v_0`$ is a finite Radon measure. Here $`v`$ is the fluid velocity, $`p`$ is the pressure, and the density of the fluid is assumed to be one. Notice that a vortex sheet gives a measure valued vorticity supported on the interface. This approach was posed by DiPerna and Majda in 1987 . In 1991, J.M. Delort  proved the existence of weak solutions global in time of the 2-D incompressible Euler equation (1) for measure-valued initial vorticity in $`H_{loc}^1(R^2)`$ that has a distinguished sign. However the problem of uniqueness of the weak solution is still unresolved. In 1963, Yudovich  obtained the existence and uniqueness of weak solutions of the 2-D incompressible Euler equation (1) for initially bounded vorticity. The best results on uniqueness upto date are given by Yudovich  and Vishik  for weak solutions with vorticity in a class slightly larger than $`L^{\mathrm{}}`$. This does not include vortex sheets, which admit measure-valued vorticity. Examples of weak solutions with the velocity field $`vL^2(R^2\times (T,T))`$ that is compactly supported in space-time was constructed by V. Scheffer  and later by A. Shnirelman . This gives non-uniqueness of weak solutions in $`L^2(R^2\times (T,T))`$. However non-uniqueness in the physically relevant class of conserved energy $`vL^{\mathrm{}}([0,\mathrm{}),L_{loc}^2(R^2))`$ remains open. Numerical evidences of non-uniqueness of weak solutions for vortex sheet data can be found in , .

Furthermore, weak solutions give little information of the specific nature of the vortex sheet evolution. For instance, does the vorticity remain supported on a curve for a later time given that the initial vorticity is supported on a curve in $`R^2`$? Assume further that the free interface between the two fluid domains remains a curve in $`R^2`$ at a later time, equation (1) can be reduced to an evolutionary differential-integral equation along the interface. This is the Birkhoff-Rott equation, written explicitly by Birkhoff in  and implied in the work of Rott . The second approach uses the Birkhoff-Rott equation as a model for the evolution of the vortex sheet.

## 2. The Birkhoff-Rott equation

For convenience, we use complex variable $`z=x+iy`$ to denote a point in $`R^2`$. $`\overline{z}=xiy`$ denotes the complex conjugate and $`f_x=_xf`$ is the partial derivative of the function $`f`$. $`H^s`$ indicates Sobolev spaces.

In search of the equation for the evolution of the vortex sheet, we suppose that at time $`t0`$ the vorticity is a measure supported on the curve $`\mathrm{\Gamma }(t)`$ given by the complex position $`\xi =\xi (s,t)`$ in the arclength $`s`$, in which $`\xi (0,t)`$ is the particle path of a reference particle; and on this curve the vorticity density is $`\gamma =\gamma (s,t)`$. That is the vorticity at time $`t`$ is $`\omega (x,y,t)`$ satisfying

$$\varphi (x,y)\omega (x,y,t)𝑑x𝑑y=\varphi (\xi (s,t))\gamma (s,t)𝑑s,\text{for any }\varphi C_0^{\mathrm{}}(R^2).$$

From the Biot-Savart law, the velocity field $`v`$ induced by the vorticity is given by

$$\overline{v}(z,t)=\frac{1}{2\pi i}\frac{\gamma (s^{},t)}{z\xi (s^{},t)}𝑑s^{},\text{for }z\mathrm{\Gamma }(t).$$

Notice that the velocity is discontinuous just on $`\mathrm{\Gamma }(t)`$. We define the velocity on the sheet as the average of the velocities at the two sides of the sheet, that is given by the principle value integral:

$$\overline{v}(\xi (s,t),t)=\frac{1}{2\pi i}p.v.\frac{\gamma (s^{},t)}{\xi (s,t)\xi (s^{},t)}𝑑s^{}.$$
$`2`$
As suggested by the properties of the Euler equation, we assume that the vortex sheet is convected by the average velocity (2), and the vorticity is conserved along the particle path. We arrive at the evolution equation of the vortex sheet:

$`\xi _t(s,t)+a(s,t)\xi _s(s,t)`$ $`=v(\xi (s,t),t)`$ $`3`$
$`\gamma _t(s,t)+_s(a(s,t)\gamma (s,t))`$ $`=0`$
where $`a(s,t)`$ is a real valued function satisfying $`a(0,t)=0`$.

A rigorous justification of the equivalence between equation (3) and equation (1) for smooth graphs $`\xi (s,t)`$ and smooth vortex strength $`\gamma (s,t)`$ can be found in . It is not hard to extend it to all smooth curves.

Assume $`\alpha (s,t)=_0^s\gamma (s^{},t)𝑑s^{}`$ defines an increasing function of $`s`$, we make a change of variables: $`z(\alpha ,t)=\xi (s(\alpha ,t),t)`$, in which $`s(\alpha ,t)`$ is the inverse of $`\alpha (s,t)`$: $`\alpha (s(\alpha ,t),t)=\alpha `$. We get from equation (3) the Birkhoff-Rott equation

$$_t\overline{z}(\alpha ,t)=\frac{1}{2\pi i}p.v.\frac{1}{z(\alpha ,t)z(\beta ,t)}𝑑\beta .$$
$`4`$
Notice that $`z=z(\alpha ,t)`$ is a parameterization of the vortex sheet in the circulation variable $`\alpha `$, $`1/|z_\alpha |=\gamma `$ is the vortex strength. A steady solution of (4) is the flat sheet $`z=\alpha `$.

Equation (4) has been under active investigations over the last four decades. A well-known property of (4) is that perturbations of the flat sheet grow due to the Kelvin-Helmholtz instability, following from a linearization of equation (4) about the flat sheet. For given analytic data, Sulem, Sulem, Bardos and Frisch  established the short time existence and uniqueness of solutions in analytic class for 2-D and 3-D vortex sheet evolution. Duchon and Robert  obtained the global existence of solutions of equation (4) for a special class of initial data that is close to the flat sheet. However, numerous results show that a vortex sheet can develop a curvature singularity in finite time from analytic data. D.W. Moore  was the first to provide analytical evidence that predicts the occurrence and time of singularity formation, which was verified numerically by Meiron, Baker and Orszag  and by Krasny . Caflisch and Orellana  proved existence almost up to the time of expected singularity formation for analytic data that is close to the flat sheet. Duchon and Robert  and Caflisch and Orellana  constructed specific examples of solutions of equation (4) where a curvature singularity develops in finite time from analytic data. The example of Caflisch and Orellana  has the form $`z(\alpha ,t)=\alpha +S(\alpha ,t)+r(\alpha ,t)`$, where

$$S(\alpha ,t)=ϵ(1i)\{(1e^{t/2i\alpha })^{1+\mu }(1e^{t/2+i\alpha })^{1+\mu }\}$$

is a solution of the linearized equation in which $`ϵ`$ is small, $`\mu >0`$; $`r(\alpha ,t)`$ is the correction term that is negligible relative to $`S(\alpha ,t)`$ in the sense that $`S(\alpha ,t)+r(\alpha ,t)`$ exhibits the same kind of behavior as $`S(\alpha ,t)`$ . Notice that $`S(\alpha ,t)`$ is an analytic function for $`t>0`$, but $`S(\alpha ,0)`$ has an infinite second derivative at $`\alpha =0`$ for $`\mu (0,1)`$. In fact, the $`(1+\nu )`$th derivative of $`S(\alpha ,0)`$ for $`\nu >\mu `$ becomes infinite at $`\alpha =0`$. Now inverting time gives an example $`\widehat{z}(\alpha ,t)`$ that is analytic at $`t_0<0`$, but has an infinite second derivative at $`\alpha =0`$, $`t=0`$. At the singularity formation time $`t=0`$, the vortex strength $`1/|\widehat{z}_\alpha |`$ of this example satisfies

$$0<c1/|\widehat{z}_\alpha |C<\mathrm{}$$
$`5`$
for some constants $`c`$ and $`C`$; and $`\widehat{z}(\alpha ,t)C^{1+\rho }(R\times [t_0,0])`$ for $`0<\rho <\mu `$.

These examples also show that the initial value problem of the Birkhoff-Rott equation (4) is ill-posed in $`C^{1+\nu }(R)`$, $`\nu >0`$, and in Sobolev spaces $`H^s(R)`$, $`s>3/2`$ in the Hadamard sense , . Ill-posedness was also proved by Ebin  using a different approach. However the existence of solutions in spaces less regular than $`C^{1+\nu }(R)`$ or $`H^s(R)`$, and the nature of the vortex sheet at and beyond the singularity time remained unknown analytically in general.

This suggests that we look for solutions of the Birkhoff-Rott equation in the largest possible spaces where the equation makes sense. For the purpose of this paper, we consider functions $`z(\alpha ,t)`$ so that for each fixed time $`t`$, both sides of the equation (4) are functions locally in $`L^2`$, and on which the $`L^2`$-analysis is available. This leads us to consider chord-arc curves, thanks to the work of G. David .

Another reason that chord-arc curves are to be considered is due to the numerical calculation of Krasny  . Krasny studied the evolution of the vortex sheet beyond singularity using the vortex blob method. He found that the approximating solutions have the form of a spiral beyond singularity. Convergence of the approximating sequence to a weak solution of the Euler equation (1) was proved by J-G. Liu and Z-P. Xin , under the assumption that the initial vorticity has a distinguished sign. A special example of chord-arc curves is a logarithmic spiral.

## 3. Chord-arc curves and some recent results

Let $`\mathrm{\Gamma }`$ be a rectifiable Jordan curve in $`R^2`$ given by $`\xi =\xi (s)`$ in the arclength $`s`$. We say $`\mathrm{\Gamma }`$ is a chord-arc curve, if there is a constant $`M1`$, such that

$$|s_1s_2|M|\xi (s_1)\xi (s_2)|,\text{for all }s_1,s_2.$$

The infimum of all such constants M is called the chord-arc constant.

For a chord-arc curve $`\xi =\xi (s)`$, $`s`$ the arclength, it is proved in  that $`\xi ^{}(s)`$ exists almost everywhere, and there is a choice of the argument function $`bBMO`$, with $`\xi ^{}(s)=e^{ib(s)}`$. In particular, if the chord-arc constant is close to 1, there is a choice of $`bBMO`$, such that $`b_{BMO}`$ is close to $`0`$. Moreover the subset of all those functions $`b`$ is an open subset of $`BMO`$. And if $`bBMO`$, and $`b_{BMO}<1`$, $`\xi (s)=\xi _0+_0^se^{ib(s^{})}𝑑s^{}`$ defines a chord-arc curve.

Examples of chord-arc curves include Lipschitz curves and logarithmic spirals $`r=\pm e^\theta `$, $`\theta R`$, where $`(r,\theta )`$ is the polar coordinates.

A Theorem of G. David  states that

###### Theorem (G. David \cite{7})

For all chord-arc curves $`\mathrm{\Gamma }:\xi =\xi (s)`$, $`s`$ the arclength, the corresponding Cauchy integral operator $`C_\mathrm{\Gamma }`$, where

$$C_\mathrm{\Gamma }f(s)=p.v.\frac{f(s^{})}{\xi (s)\xi (s^{})}𝑑\xi (s^{}),$$

is bounded from $`L^2(ds)`$ to $`L^2(ds)`$.

In fact, the result of G. David  is stronger than stated above. He proved that the Cauchy integral operator $`C_\mathrm{\Gamma }`$ is bounded from $`L^2(ds)`$ to $`L^2(ds)`$ if and only if $`\mathrm{\Gamma }`$ is a regular curve. A rectifiable curve $`\mathrm{\Gamma }`$ is said to be regular if there is a constant $`M`$ such that for every $`r>0`$ and every disc $`D`$ with radius $`r`$, the length of $`\mathrm{\Gamma }D`$ does not exceed $`Mr`$. A chord-arc curve is regular but not vice versa.

Now we go back to the Birkhoff-Rott equation. Notice that the Biot-Savart integral representing an incompressible velocity field $`v`$ in terms of the vorticity $`\omega `$ may be divergent if $`\omega `$ does not vanish fast enough at infinity, even if the velocity field $`v`$ is well defined, we extend the definition of the Birkhoff-Rott equation (4) by considering the differences of the velocities between any two points:

$`\overline{z_t}(\alpha ,t)`$ $`\overline{z_t}(\alpha ^{},t)={\displaystyle \frac{1}{2\pi i}}p.v.{\displaystyle _{|\beta |N}}\{{\displaystyle \frac{1}{z(\alpha ,t)z(\beta ,t)}}{\displaystyle \frac{1}{z(\alpha ^{},t)z(\beta ,t)}}\}𝑑\beta `$ $`6`$
$`+{\displaystyle \frac{1}{2\pi i}}p.v.{\displaystyle _{|\beta |>N}}{\displaystyle \frac{z(\alpha ^{},t)z(\alpha ,t)}{(z(\alpha ,t)z(\beta ,t))(z(\alpha ^{},t)z(\beta ,t))}}𝑑\beta ,`$
for all $`(\alpha ,t)`$, $`(\alpha ^{},t)`$, and some $`N>|\alpha |+|\alpha ^{}|+1`$. This admits a larger class of solutions. In particular, the integral on the right hand side of equation (6) is convergent for those similarity solutions considered in , -, which otherwise give divergent Cauchy integrals in (4) due to the divergent contributions from the vorticities at infinity. It follows from the Theorem of G. David that the integral on the right hand side of the equation (6) is convergent for $`a.e.`$ $`(\alpha ,t)`$, $`(\alpha ^{},t)`$ and is in $`L^{\mathrm{}}([0,T],L_{loc}^2(d\alpha )\times L_{loc}^2(d\alpha ^{}))`$ for the solutions considered in Theorem 1 in the following.

Roughly speaking, if a function $`z=z(\alpha ,t)`$ satisfies equation (4), it will also satisfy equation (6). On the other hand, if $`z=z(\alpha ,t)`$ satisfies (6), and the Cauchy integral

$$\frac{1}{2\pi i}p.v.\frac{1}{z(\alpha ,t)z(\beta ,t)}𝑑\beta $$

is convergent, then there is a function $`c=c(t)`$, such that $`z(\alpha ,t)+c(t)`$ satisfies equation (4).

For a local integrable function $`f=f(\alpha )`$ defined on $`(a,b)`$, we say $`f`$ is of bounded local mean oscillation on $`(a,b)`$ if there exists $`\delta _0>0`$ such that

$$f_{BMO(a,b),\delta _0}=\underset{\text{all }I(a,b),|I|\delta _0}{sup}\frac{1}{|I|}_I|f(\alpha )f_I|𝑑\alpha <\mathrm{},$$

here $`f_I=\frac{1}{|I|}_If(\alpha )𝑑\alpha `$, $`I`$ is an interval. We say $`f`$ is analytic on $`(a,b)`$ if $`fC^{\mathrm{}}(a,b)`$, and for any compact subset $`K`$ of $`(a,b)`$, there is a constant $`\rho >0`$, such that

$$\underset{m=o}{\overset{\mathrm{}}{}}\frac{\rho ^m}{m!}_K|_\alpha ^mf(\alpha )|^2𝑑\alpha <\mathrm{}.$$

Notice that $`\mathrm{ln}z`$ is multi-valued for complex number $`z`$. In the following, $`\mathrm{ln}z_\alpha `$ refers to one choice of the multi-values. We have the following results concerning the solutions of the Birkhoff-Rott equation (6) (or (4)).

###### Theorem 1 \cite{32}

Assume that $`zH^1([0,T],L_{loc}^2(R))L^2([0,T],H_{loc}^1(R))`$ is a solution of the Birkhoff-Rott equation (6) for $`0tT`$, satisfying that

1. There are constants $`m>0`$, $`M>0`$, independent of $`t`$, such that

$$m|\alpha \beta ||z(\alpha ,t)z(\beta ,t)|M|\alpha \beta |\text{for all }\alpha ,\beta ,0tT.$$
$`7`$
2. For the interval $`(a,b)`$, there exists $`\delta _0>0`$, independent of $`t`$, such that for all $`0tT`$,

$$\mathrm{ln}z_\alpha (,t)_{BMO(a,b),\delta _0}C(m,M),$$
$`8`$
here $`C(m,M)`$ is a universal constant depending on $`m`$ and $`M`$. Assume further that $`\mathrm{ln}z_\alpha L^2([0,T],L_{loc}^2(R))`$. Then $`z_\alpha C((a,b)\times (0,T))`$, and for each $`t_0(0,T)`$, $`z_\alpha (,t_0)`$ is analytic on $`(a,b)`$.

Remark  Notice that assumption 2. is satisfied if $`\mathrm{ln}z_\alpha C([a,b]\times [0,T])`$. Assumption 1. is equivalent to assuming that the vortex strength $`\gamma =1/|z_\alpha |`$ is bounded away from $`0`$ and $`\mathrm{}`$, with bounds independent of $`t`$, and $`z=z(,t)`$ defines a chord-arc curve for each fixed $`t[0,T]`$, with the chord-arc constant independent of $`t`$. Assumption 1. can be relaxed by requiring that $`1/|z_\alpha (\alpha ,t)|`$ is bounded away from $`0`$ and $`\mathrm{}`$ for $`\alpha (a,b)`$ only, and by assuming some weaker conditions for $`\alpha (a,b)`$.

A version of Theorem 1 under assumptions that the solution of the Birkhoff-Rott equation (4): $`z=z(\alpha ,t)C^{1+\rho _0}`$, for some $`\rho _0>0`$, the vortex strength $`c_0<\gamma (\alpha ,t)<C_0`$ for some constants $`c_0>0`$, $`C_0>0`$, and the vortex sheets $`\mathrm{\Gamma }(t)`$ are closed Jordan curves is also obtained by Lebeau  using an independent approach.

As a consequence of Theorem 1, the example constructed by Caflisch and Orellana  will fail to satisfy properties 1 and 2 as stated in Theorem 1 near the singularity after the singularity formation time.

Let $`\text{Re}z`$ and $`\text{Im}z`$ be the real and imaginary parts of the complex number $`z`$ respectively. Regarding the existence of solutions, we have the following

###### Theorem 2 \cite{32}

For any real valued function $`w_0H^{\frac{3}{2}}(R)`$, there exists $`T=T(w_0_{H^{\frac{3}{2}}})>0`$, such that the Birkhoff-Rott equation (6) has a solution $`z=z(\alpha ,t)`$ for $`0tT`$, satisfying $`\mathrm{ln}z_\alpha L^{\mathrm{}}([0,T],H^{\frac{3}{2}}(R))Lip([0,T],H^{\frac{1}{2}}(R))`$ and $`\text{Im}\{(1+i)\mathrm{ln}z_\alpha (\alpha ,0)\}=w_0(\alpha )`$, with the properties that there exist constants $`m>0`$, $`M>0`$ independent of $`t`$, such that

$$m|\alpha \beta ||z(\alpha ,t)z(\beta ,t)|M|\alpha \beta |,\text{for all }\alpha ,\beta ,0tT;$$

and there exists $`\delta _0>0`$, independent of $`t`$, such that

$$\mathrm{ln}z_\alpha (,t)_{BMO(R),\delta _0}C(m,M),\text{for }0tT,$$

here $`C(m,M)`$ is the universal constant as in Theorem 1.

Theorem 2 states that if only half of the data $`z_\alpha (\alpha ,0)`$ is given, there is a solution of the Birkhoff-Rott equation (6) for a finite time period. Theorem 2 is a generalization of the existence result of Duchon and Robert  to general data.

The following result implies that Theorem 2 is optimal, in the sense that in general, there is no solution of the Birkhoff-Rott equation satisfying properties 1. and 2. as stated in Theorem 1 beyond the initial time $`t=0`$ for arbitrarily given data.

###### Theorem 3 \cite{32}

Assume that $`zH^1([0,T],L_{loc}^2(R))L^2([0,T],H_{loc}^1(R))`$ is a solution of the Birkhoff-Rott equation (6) for $`0tT`$, $`T>0`$, satisfying the property 1. and property 2. on some interval $`(a,b)`$ as stated in Theorem 1. Assume further that $`\mathrm{ln}z_\alpha L^2([0,T],L_{loc}^2(R))`$, and $`w_0=\text{Im}\{(1+i)\mathrm{ln}z_\alpha (,0)\}`$ is analytic on $`(a,b)`$. Then $`z_\alpha C((a,b)\times [0,T))`$ and $`\text{Re}\{(1+i)\mathrm{ln}z_\alpha (,0)\}`$ is also analytic on $`(a,b)`$.

## 4. Open questions

The reason that Theorem 1-3 holds is that under the assumptions 1. and 2. in Theorem 1, the Birkhoff-Rott equation (6) is of “elliptic” type on $`(a,b)\times [0,T]`$. It would be interesting to see whether the assumption 2 that requires small local mean oscillation on $`\mathrm{ln}z_\alpha (,t)`$ can be removed, or to see whether one can construct a non-smooth solution of equation (6) that violates (8). A good place to start is to construct similarity solutions of the Birkhoff-Rott equation. Similarity solutions are studied numerically in Pullen  , Pullen and Phillips , and analytically in Kambe . However, the similarity solutions in , - violate (7). It would also be interesting to relax the assumption 1. by considering vortex strength that is not necessarily bounded away from $`0`$ and infinity on the interval $`(a,b)`$.

Theorem 1-3 implies that a solution of the Birkhoff-Rott equation is either analytic, or doesn’t satisfy properties 1. and 2. in Theorem 1. And there are in general no solutions of reasonable regularity (properties 1. and 2.) beyond the initial time for an arbitrarily given data. Are there solutions of the Birkhoff-Rott equation in spaces of even less regularity? Notice that we can define the solutions of the Birkhoff-Rott equation in the distribution sense by

$$_t(\overline{z}(\alpha ,t)\eta (\alpha )𝑑\alpha )=\frac{1}{4\pi i}\frac{\eta (\alpha )\eta (\beta )}{z(\alpha ,t)z(\beta ,t)}𝑑\alpha 𝑑\beta ,\text{for all }\eta C_0^{\mathrm{}}.$$

This admits some even larger classes of solutions $`z=z(\alpha ,t)`$. Is there a solution in the distribution sense for the Birkhoff-Rott equation for any given data? How is it related to the weak solution of the Euler equation?

Theorem 1-3 might as well suggest that the vortex sheet in general fails to be a curve beyond the initial time for general data. Therefore it becomes interesting to study the vortex layers or considering the effects of viscosity. Numerical analysis of the vortex layers can be found in Baker and Shelley  etc.

## References


format: text

text:
# Zeeman levels with exotic field dependence in the high field phase of a 𝑺=1 Heisenberg antiferromagnetic chain

## Abstract

We have performed electron spin resonance measurements over a wide frequency and magnetic field range on a single crystal of the $`S=1`$ quasi-one-dimensional Heisenberg antiferromagnet $`\mathrm{Ni}(\mathrm{C}_5\mathrm{H}_{14}\mathrm{N}_2)_2\mathrm{N}_3(\mathrm{PF}_6)`$. We observed gapped excitation branches above the critical field ($`H_\mathrm{c}`$) where the Haldane gap closes. These branches are analyzed by a phenomenological field-theory using complex-field $`\varphi ^4`$ model. A satisfactory agreement between experiment and theory is obtained.

Recently, field-induced phenomena in quantum spin systems have attracted considerable interest. These include magnetization plateaus narumi ; shiramura ; kageyama  and field-induced long-range order (LRO) honda1 ; honda2 ; chen . Gapped one-dimensional (1D) spin systems subject to an external magnetic field strong enough to close the gap are driven into a new phase. When the system has XY or Heisenberg symmetry, this phase is critical and its low-energy physics is described by the Tomonaga-Luttinger liquid chitra , whose elementary excitations are of the particle-hole type (spinon pairs). The high-energy physics which cannot be described using this picture has been recently investigated experimentally orendac ; zheludev ; ruegg  as well as theoretically KolezhukMikeska .

One-dimensional Heisenberg antiferromagnets (HAFs) with integer spin are typically gapped and have a disordered spin liquid singlet ground state haldane . The first excited state is a triplet and one of the Zeeman-split triplet branches goes down in magnetic field. At some critical field $`H_\mathrm{c}`$ the gap closes, a qualitatively new ground state emerges and magnetism recovers. Above $`H_\mathrm{c}`$, a LRO is expected to occur in quasi-1D integer-spin HAFs owing to interchain interaction. Such a field-induced LRO in the $`S=1`$ quasi-1D HAF was found recently by specific heat measurements on a single crystal of $`\mathrm{Ni}(\mathrm{C}_5\mathrm{H}_{14}\mathrm{N}_2)_2\mathrm{N}_3(\mathrm{PF}_6)`$, abbreviated NDMAP honda1 ; honda2 . Figure 1 shows the magnetic field versus temperature phase diagram of NDMAP. The regions below and above the boundary denoted as $`H_{\mathrm{LRO}}(T)`$ correspond to the Haldane disordered phase, and to the LRO phase, respectively. The critical field $`H_\mathrm{c}`$ in the present paper is defined as $`H_{\mathrm{LRO}}`$ extrapolated to $`T=0`$.

Quite recently, three distinct excitation branches were observed in inelastic neutron scattering (INS) experiments on NDMAP in the LRO phase above $`H_\mathrm{c}`$ for $`Ha`$, and a satisfactory desription within a phenomenological field theory was obtained zheludev . This observation is remarkable because one would expect two branches from the conventional spin-wave theory.

Electron spin resonance (ESR) is one of the powerful methods to investigate magnetic excitations, especially at low energies with a much better energy resolution than INS. In Ref. honda3, , we observed ESR signals in a single crystal of NDMAP and interpreted the results as the coexistence of 1D and 3D excitations, because no satisfactory description was available above $`H_\mathrm{c}`$ at that time.

In the present paper we extend the ESR study of NDMAP to a much wider frequency range from which we obtain new results. A nontrivial field dependence of the resonance modes, very different from that expected from a conventional spin-wave theory, is observed. We compare the experimental results with model calculations based on the field theory used in the analysis of the recent INS experiments zheludev , and obtain a good agreement with the theory.

First we summarize the crystal and magnetic properties of NDMAP. This compound crystallizes in the orthorhombic system (space group $`Pnmn`$monfort . The $`\mathrm{Ni}^{2+}`$ ions are bridged by azido ions forming antiferromagnetic chains along the $`c`$ axis. All the $`\mathrm{Ni}`$ sites in a chain are crystallographically equivalent, so that no staggered components of the magnetic moments are retained in the ground state below $`H_\mathrm{c}`$. The physics in NDMAP is, therefore, essentially different from that in another typical Haldane chain material $`\mathrm{Ni}(\mathrm{C}_2\mathrm{H}_8\mathrm{N}_2)_2\mathrm{NO}_2(\mathrm{ClO}_4)`$, abbreviated NENP, extensively studied by ESR lu-sieling-brill . These studies showed that the Haldane gap did not close at Hc due to a staggering of the local crystal fields for Ni2+MitraHalperin94 ; SakaiShiba94 . Consequently, in NENP no clear transition from the disordered to LRO phase is observed and the ESR transitions from the ground state are always allowed.

From the analysis of the magnetic susceptibility data, the following values are obtained honda1 : $`J/k_\mathrm{B}=30.0`$ K, $`D/J=0.3`$, $`g_{}=2.10`$ and $`g_{}=2.17`$, where $`J`$ is the intrachain exchange interaction constant and $`g_{}`$ and $`g_{}`$ are the $`g`$ values parallel and perpendicular to the chain, respectively. From the INS measurements done on single crystals of deuterated NDMAP zheludev2  the following parameters were determined: $`J=2.85`$ meV($`=33.1`$ K), $`J_x^{}=3.5\times 10^4`$ meV ($`=4.1\times 10^3`$ K), $`J_y^{}=1.8\times 10^3`$ meV ($`=2.1\times 10^2`$ K) and $`D=0.70`$ meV ($`=8.1`$ K), where $`J_x^{}`$ and $`J_y^{}`$ are the interchain exchange interactions along the $`a`$ and $`b`$ axes, respectively.

The single crystals of NDMAP were grown in the same manner as that described in Ref.honda1, . The ESR measurements were performed using a millimeter vector network analyzer with extensions up to about 700 GHz from the AB Millimeter, France and a superconducting magnet up to 16 T from the Oxford Instruments, UK. We performed the ESR measurements in $`Ha`$ and $`Hc`$ geometries at the lowest temperature available with our spectrometer ($`1.5`$ K). Faraday and Voigt configurations were used for $`Ha`$ and $`Hc`$, respectively. As shown in Fig. 1, the external magnetic field of $`9`$ T for $`Ha`$ brings the system into the LRO phase at $`1.5`$ K.

Figure 2 shows a typical example of ESR spectra observed for $`Ha`$ and $`Hc`$. The signal at about $`13.5`$ T for $`Ha`$ comes from the LRO phase. In Fig. 3(a), the resonance frequencies for $`Ha`$ as a function of magnetic field are shown with solid circles. We also plot with solid squares the gap energies obtained from the recent INS measurements at 30 mK zheludev . In the figure, the vertical broken line and the vertical dotted line denote the $`T=0`$ critical field $`H_c`$ and the boundary field $`H_{\mathrm{LRO}}`$ between the disordered and LRO phases at $`1.5`$ K, respectively. It is obvious that the slope of the branches changes at $`H_\mathrm{c}`$ and not at $`H_{\mathrm{LRO}}`$. It should be remarked that a similar change of slope was observed in the quasi-1D material NENC orendac  and in the 3D weakly coupled dimer system $`\mathrm{TlCuCl}_3`$ ruegg ; Matsumoto+02 . The observed neutron data below $`H_\mathrm{c}`$ represent the magnon triplet split by the single ion anisotropy in the Haldane phase. Below $`H_\mathrm{c}`$, only temperature-induced ESR transitions between the magnon branches are allowed. Only one ESR line, corresponding to the transition between the lowest two branches, was observed near the critical field below $`H_\mathrm{c}`$. Above $`H_\mathrm{c}`$, three ESR lines were observed. Their positions coincide with the INS data for the gaps at $`q=\pi `$, confirming that they correspond to the transitions from the ground state. It is worthwhile to note that these lines appear at $`H>H_\mathrm{c}`$ and not at $`H>H_{\mathrm{LRO}}`$ as one would naively expect, and that there is no anomaly in the behavior of resonance lines at $`H_{\mathrm{LRO}}`$. Indeed, the gap of the lowest magnon branch at $`q=\pi `$ closes at $`H=H_\mathrm{c}`$ at any temperature, causing “condensation” of $`q=\pi `$ magnons into the ground state at $`H>H_\mathrm{c}`$. The amplitude of the $`q=\pi `$ component becomes nonzero immediately above $`H_\mathrm{c}`$, making the ESR transitions from the ground state to the excited ones at $`q=\pi `$ possible. The transversal staggered LRO, however, emerges only at the higher field $`H_{\mathrm{LRO}}`$, because the phase coherence of the condensate wave function, directly related to the LRO, is destroyed by thermal fluctuations in the field range $`H_\mathrm{c}<H<H_{\mathrm{LRO}}`$. When the phase coherence sets in at $`H=H_{\mathrm{LRO}}`$, there is no peculiarity in the amplitude and thus no essential change in the ESR signal.

Figure 3(b) shows a frequency versus field plot for $`Hc`$. As in the $`Ha`$ case, we observed only one resonance line below $`H_\mathrm{c}`$ which corresponds to the temperature-induced transitions between the lowest two magnon branches. Above $`H_\mathrm{c}`$, two ESR lines were observed. For this geometry, there is no LRO at $`1.5`$ K, so that the measurements were conducted completely outside the LRO phase. One of the remarkable features in this figure is a change in slope at about 7 T, far above the critical field $`H_\mathrm{c}3.5`$ T, resulting from an avoided crossing between the two branches (see below).

Now, we compare the experimental results with theoretical calculations using a field-theoretical model zheludev  similar to that proposed for dimerized $`S=\frac{1}{2}`$ chains in Ref. K96, . It is a Ginzburg-Landau-type theory written in terms of a *complex* triplet field $`\mathit{\varphi }=𝑨+i𝑩`$ which is assumed to be small, $`|\mathit{\varphi }|1`$. For a dimerized $`S=\frac{1}{2}`$ system, the field $`\mathit{\varphi }`$ describes the dimer wave function $`|\mathit{\varphi }=(1|\mathit{\varphi }|^2)^{1/2}|s+_{a=x,y,z}\varphi _a|t_a`$, where $`|s`$ and $`|t_a`$ are the singlet and three triplet states of the dimer; it is a continuum analog of the bond boson description used in Matsumoto+02 . The uniform and staggered magnetization are given by $`𝑴=(\mathit{\varphi }^{}\times \mathit{\varphi })`$ and $`𝑳=(\mathit{\varphi }^{}+\mathit{\varphi })(1|\mathit{\varphi }|^2)`$, respectively. The effective Lagrangian density in the continuum limit has the following form:

$``$ $`=`$ $`\mathrm{}(𝑨_t𝑩𝑩_t𝑨)(1/2)v^2(_x𝑨)^2`$
$``$ $`{\displaystyle \underset{\alpha }{}}\{m_\alpha A_\alpha ^2+\stackrel{~}{m}_\alpha B_\alpha ^2\}+2𝑯(𝑨\times 𝑩)`$
$``$ $`\lambda (𝑨^2)^2\lambda _1(𝑨^2𝑩^2)\lambda _2(𝑨𝑩)^2.`$
This model, though derived for dimerized $`S=\frac{1}{2}`$ chains, may also be used for other gapped 1D systems which are in the same universality class, particularly for $`S=1`$ Haldane chains. In this case, the quantities $`m_\alpha `$, $`\stackrel{~}{m_\alpha }`$, and $`\lambda _i`$ should be treated as phenomenological parameters; the direct microscopic justification of this model for $`S=1`$ chains is the subject of ongoing work. The spatial derivatives of $`𝑩`$ are omitted in (Zeeman levels with exotic field dependence in the high field phase of a $`𝑺=1`$ Heisenberg antiferromagnetic chain) because they appear only in terms which are of the fourth order in $`𝑨`$ and $`𝑩`$. By integrating out the “slave” $`𝑩`$-field one obtains an effective real-field $`\phi ^4`$-type theory similar to that of Affleck Affleck90-91 ; actually, Affleck’s theory may be viewed as a special case of (Zeeman levels with exotic field dependence in the high field phase of a $`𝑺=1`$ Heisenberg antiferromagnetic chain) with isotropic $`\stackrel{~}{m}_\alpha =\stackrel{~}{m}`$ and with simplified interaction term $`\lambda _{1,2}=0`$. Another special case $`m_\alpha =\stackrel{~}{m_\alpha }`$ leads to the spectra which for $`H<H_c`$ exactly coincide with those obtained in the approach of Tsvelik Tsvelik90  who has proposed an effective theory involving three massive Majorana fields, and also with the perturbative formulas of Refs. perturbative, . If one additionally assumes the simplified interaction potential $`\lambda _{1,2}=0`$, the theory becomes equivalent to that of Mitra and Halperin MitraHalperin94  who postulated a bosonic Lagrangian to match Tsvelik’s results for the field dependence of the gaps below $`H_c`$.

The solid lines in Figs. 3(a) and (b) represent the best fit of the calculated energy gaps at $`q=\pi `$ to the experimental data. This fit is obtained using the following set of model parameters (all values in meV): $`m_a=0.5`$, $`m_b=0.8`$, $`m_c=5.29`$, $`\stackrel{~}{m_a}=0.353`$, $`\stackrel{~}{m_b}=0.342`$, $`\stackrel{~}{m_c}=0.647`$, $`\lambda =1`$, $`\lambda _1=\lambda _2=0.13`$. We obviously see a change in slope of each excited branch and re-opening the energy gap of the lowest branch at $`H_\mathrm{c}`$. Below $`H_\mathrm{c}5.5`$ T in Fig. 3(a), three excitation branches are in a good agreement with the neutron data. The ESR transitions from the ground to the excited states at $`q=\pi `$ are forbidden because no staggered component exists in the ground state of this compound. In both figures, the ESR resonance points only appear near $`H_\mathrm{c}`$, because these correspond to the transition between the excited states and the intensity is related to the population of the excited states at $`1.5`$ K. The thin broken line in each figure indicates the transition between the lowest two magnon branches at $`q=\pi `$. For those temperature-induced transitions, the agreement between experiment and calculation is satisfactory for $`Ha`$ and is not good for $`Hc`$. Above $`H_\mathrm{c}`$, the calculated branches are close to the experimental data near $`H_\mathrm{c}`$, but deviate from them at high fields. This is probably because the theory is applicable only for the state with a small staggered order and thus only for fields which are not too far above $`H_\mathrm{c}`$.

In these calculations, we have taken into account the tilting of the crystal field axis of the $`\mathrm{Ni}^{2+}`$ ions from the crystallographic $`c`$ axis. The tilting angle is $`15.9^{}`$. In Fig. 3(b), a change in slope at about 7 T observed in the A and B branches are reproduced in the calculation. Because of the tilting of the crystal field axis, the external magnetic field applied parallel to the $`c`$ axis mixes the wave functions of A and B branches. The origin of this change in slope is different from that of the changes visible near the critical field. The observed ESR line below $`H_\mathrm{c}`$ in Fig. 3(b) deviates from the expected thin broken line, although a satisfactory agreement is obtained for $`Ha`$; the reason for this discrepancy is not clear at present. For $`Hc`$, we did not observe ESR signals corresponding to the lowest calculated mode above $`H_\mathrm{c}`$; we believe this is caused by a strong damping of the lowest magnon mode due to the interaction with domain walls, because our measurement was done at a high temperature. A similar problem arose in INS measurements Zheludev+02-PRL  in the $`Ha`$ geometry done at $`2`$ K where the lowest magnon mode appeared to be quasielastic, and only after lowering the temperature to $`30`$ mK the mode was successfully observed zheludev . The damping should be much stronger in the $`Hc`$ geometry, due to the lower energy and respectively higher density of domain walls.

Remarkably, the observed field dependence of the ESR modes for both geometries is very different from that expected from a conventional spin-wave theory. Indeed, in a conventional (classical) antiferromagnet with easy-plane anisotropy $`D`$ and external magnetic field $`H`$ applied in the easy plane, there are two resonance modes with the energies $`\epsilon _1=g\mu _BH`$ and $`\epsilon _2=2S\sqrt{2DJ}`$. For $`H`$ perpendicular to the easy plane, one would classically expect to see one field-dependent resonance mode with the energy $`\epsilon _2^{}=\sqrt{(g\mu _BH)^2+\epsilon _2^2}`$, the energy of the other mode being zero. We have observed and described theoretically completely different behaviors above $`H_\mathrm{c}`$ for both directions: in the $`Ha`$ case (in-plane geometry) the energies of all modes increase approximately as $`g\mu _BH`$, and for $`Hc`$ we observe the (avoided) crossing of two modes whose energies behave roughly as $`g\mu _BH`$ and $`2g\mu _BH`$.

In conclusion, we have performed ESR measurements on a single crystal of the $`S`$=1 quasi-one-dimensional Heisenberg antiferromagnet NDMAP and observed ESR spectra corresponding to the gapped state above the critical field where the Haldane gap closes. For $`Ha`$ our ESR data above $`H_\mathrm{c}`$, taken at $`T=1.5`$ K, agree very well with those from inelastic neutron scattering zheludev  obtained at much lower temperature, and the observed resonances are not affected by the boundary of the long-range-ordered phase at $`H9`$ T. For $`Hc`$ we observed a change in slope of the ESR line resulting from an avoided crossing between the two modes having nontrivial slopes $`H`$ and $`2H`$. The observed field dependence of the magnon modes for both geometries $`Ha`$ and $`Hc`$ is very different from that expected from a conventional spin-wave theory. Model calculations based on the phenomenological field-theoretical approach have reproduced all the features of the observation and a satisfactory quantitative agreement between theory and experiment is obtained.

*Acknowledgments.—* We would like to thank A. Zheludev for enlightening discussions. This work was in part supported by the Molecular Ensemble research program from RIKEN and the Grant-in-Aid for Scientific Research on Priority Areas(B): Field-Induced New Quantum Phenomena in Magnetic Systems (No.13130203) from the Japanese Ministry of Education, Culture, Sports, Science and Technology, and by the Grant I/75895 “Low-dimensional magnets in high magnetic fields” from Volkswagen-Stiftung.


text:
# Dark-in-Bright Solitons in Bose-Einstein Condensates with Attractive Interactions

## I Introduction

The current experimental and theoretical studies of Bose-Einstein condensates (BECs)  have attracted a great deal of interest to nonlinear patterns which can exist in them, including dark  and bright  solitons. Two-dimensional (2D) excitations in the form of vortices were also realized experimentally . Other nonlinear excitations, such as Faraday waves , ring dark solitons and vortex necklaces  were predicted to occur in BECs.

The behavior of the condensate crucially depends on the sign of the atomic interactions: dark (bright) solitons can be created in BECs with repulsive (attractive) interactions, resulting from the positive (negative) scattering length. Recent experiments have demonstrated that “tuning” of the effective scattering length, including a possibility to change its sign, can be achieved by means of a Feshbach resonance .

In this work we demonstrate that another type of localized nonlinear excitations, which may be regarded as dark solitons embedded in bright ones, can be created in attractive BECs. They may also be considered as bound states of two bright solitons with a phase shift $`\pi `$ between them. In some respects, we follow the path of very recent results of Ref. . These, in turn, were based on earlier works on the so-called twisted localized modes (TLMs) in discrete systems ; that is why we apply the same term, TLM, to these solitons of the combined dark-bright type. TLMs in a discrete system are bright solitons in which the field vanishes, changing its sign, at the central point. An alternative way to view a TLM is as a concatenation of an “up-pulse” (i.e., $`\mathrm{sech}(xx_1)`$) with a “down pulse” (i.e., $`\mathrm{sech}(xx_2)`$), where $`x_1`$ and $`x_2`$ are appropriately separated, to form an up-down, two-pulse configuration. Notice, however, that our work is essentially different from Ref. , since we focus on the case of BECs with negative scattering length (such as 7Li  and 85Rb ), and demonstrate robustness and stability of TLMs in the latter context. These structures, if regarded as effectively dark solitons in the attractive BEC, are the inverse of recently predicted  bright solitons in repulsive BECs with an optical lattice. It is expected that a characteristic longitudinal length of the TLM soliton in physical units will be on the order of the wavelength of light which induces the OL, i.e., $`1`$ $`\mu `$m, and the number of atoms in the soliton may typically be $`10^3`$$`10^4`$.

We consider a quasi-1D (cigar-shaped) BEC with negative scattering length, the corresponding normalized Gross-Pitaevskii (GP) equation being

$$iu_t=u_{xx}|u|^2u+\left[ϵx^2+V_0\mathrm{cos}^2\left(2\pi x/\lambda \right)\right]u,$$
(1)
where $`u(x,t)`$ is the macroscopic wave function, $`ϵ`$ measures the strength of the parabolic magnetic trap, and the last term accounts for the OL potential created by interference of two counterpropagating coherent light beams . As we consider only (effectively) 1D configurations, the attractive interactions between the atoms (characterized by the negative scattering length) can not produce collapse in the framework of the present model, hence various solitons have a chance to be stable.

It is obvious that scaling invariance of Eq. (1) makes it possible to fix the optical wavelength, setting $`k2\pi /\lambda =1`$. This essentially means that all the remaining length scales of the problem are measured by comparison to the period of the optical lattice which has been fixed to $`\pi `$. We follow this convention below. The number of atoms in the condensate, given by

$$N=_{\mathrm{}}^+\mathrm{}|u(x)|^2𝑑x,$$
(2)
is a dynamical invariant of Eq. (1).

In the next section we report detailed numerical results for TLM solutions in attractive BECs, as modeled by Eq. (1). In section III, we will support some of the numerical findings by analytical considerations. Section IV concludes the paper.

## II Numerical results

### A Stationary twisted-localized-mode solutions and their linear stability

We seek stationary solutions to Eq. (1) in the form,

$`u(x,t)=\mathrm{exp}(i\mu t)v(x),`$

where $`\mu `$ is the chemical potential, and a real function $`v(x)`$ obeys the equation

$$\mu v=v_{xx}v^3+\left[ϵx^2+V_0\mathrm{cos}^2(2\pi x/\lambda )\right]v\text{,}$$
(3)
supplemented with free boundary conditions. We have checked that the boundary conditions do not significantly affect the results. The corresponding boundary-value problem was solved by means of a finite-difference discretization. A second order difference scheme with spacing $`\mathrm{\Delta }x=0.2`$ was used to approximate $`v_{xx}`$. The resulting set of nonlinear algebraic equations was solved by means of a Newton iteration. More details on finite-difference schemes and Newton type methods can be found in . The solutions of interest were obtained by choosing, as the initial condition for the Newton iterations, the following functional form:

$$v(x)=\underset{j=1}{\overset{n}{}}(1)^j\eta \mathrm{sech}\left(\frac{\eta }{\sqrt{2}}(x\xi _j)\right),\eta \sqrt{\mu }.$$
(4)
The ansatz (4) approximates a superposition of nonlinear-Schrödinger (NLS) solitons. The presence of the OL suggests choosing positions of the solitons’ centers, $`\xi _j`$, as $`\xi _j=j\pi /2`$, with some integer $`j`$. The sign factor $`(1)^j`$ implies the phase difference $`\pi `$ between adjacent solitons. This choice of the phase pattern is suggested by the known result that, in the discrete NLS equation, only such a configuration may be stable, see Ref.  and references therein. Similar results for a BEC model with elliptic-function potentials were obtained in Ref. . The number $`n`$ of solitons in the initial ansatz (4) was taken as $`n=2`$, in order to find the TLM soliton proper, or $`n=3`$, with the objective to find a bound state of two TLMs, see below. However, it will become clear below that concatenating more of such “building block” structures (such as the ones with $`n=2`$ or $`n=3`$) one can, in principle, construct multi-pulse configurations with a larger number ($`n`$) of elementary pulses.

Once TLM states were found as solutions to Eq. (3) \[starting with the initial approximation (4)\], their stability was subsequently analyzed by examining the corresponding linearized problem for small perturbations. To this end, the perturbed solution was taken as

$$u(x,t)=e^{i\mu t}\left\{v(x)+\delta \left[a(x)e^{i\omega t}+b(x)e^{i\omega ^{}t}\right]\right\},$$
(5)
where $`\delta `$ is an infinitesimal amplitude of the perturbation, $`a(x)`$ and $`b(x)`$ are eigenfunctions of a perturbation mode, and $`\omega `$ is the corresponding eigenfrequency; $``$ standing for the complex conjugation. Using the ansatz of Eq. (5) to O$`(\delta )`$, one obtains the linearization problem as an eigenvalue problem where $`\omega `$ is the corresponding eigenfrequency, while $`(a,b^{})^T`$ (where the $`T`$ denotes transpose), is the corresponding eigenfunction. Using the same finite difference scheme as the one used for computing $`v(x)`$, we convert this to a matrix eigenvalue problem and find all the eigenvalues of the corresponding matrix. This is done by standard matrix eigenvalue solvers built into Matlab . $`\omega `$ can, in general, be complex, i.e., $`\omega =\omega _r+i\omega _i`$, where the subscripts denote the real and imaginary part of the eigenfrequency. If eigenfrequencies with non-zero imaginary parts exist, they lead to exponential instabilities, while their absence implies linear stability of the solution under consideration. This is monitored by the spectral plane pictures of the imaginary ($`\omega _i`$) versus the real ($`\omega _r`$) part of the eigenfrequency (see e.g., the bottom right panels of Figs. 1-3 below). The presence of even a single $`\omega `$ with $`\omega _i0`$ implies instability. In all the cases when stable stationary solitons were identified according to this criterion (see below), the stability was also checked and confirmed by direct simulations of the full equation (1). For more details on the setup and solution of the eigenvalue problem, the interested reader is referred to .

To systematically trace the evolution of numerically exact TLM solutions of Eq. (1), we performed one-parameter continuations in the three-dimensional parameter space of Eq. (3), $`(\mu ,ϵ,V_0)`$. Recall that we have set $`\lambda 2\pi `$, hence $`\lambda `$ is not a free parameter. We start by varying the OL-potential’s strength $`V_0`$ for fixed $`ϵ=0`$ and $`\mu =1.5`$. We note that negative values of the chemical potential are typical of cases when the corresponding wave-function pattern is localized, although $`\mu `$ can sometimes be positive (see below), which is explained by the presence of the OL potential in Eq. (3).

Numerical results corresponding to this case (varying $`V_0`$) are presented in Fig. 1. In the top subplot of the upper part of the figure, the number of atoms in the soliton, defined as per Eq. (2), is shown versus $`V_0`$, and the largest instability growth rate, found from the associated linearized problem, is shown in the bottom panel of the upper part of Fig. 1 (zero instability growth rate implies that the soliton is stable). Three characteristic examples of the TLM solution are displayed, together with their linear-stability eigenfrequencies, in the lower part of Fig. 1.

This solution branch exists only for $`V_0>0.075`$. The presence of this cutoff point is expected, as the optical lattice is crucial for the existence of TLM solutions, which are not found in the usual NLS equation (at least, in the absence of the parabolic-trap potential). Nevertheless, below we will present a case in which a TLM/dark soliton solution is possible even in the absence of OL, being supported solely by the magnetic trap. In the interval $`0.275<V_0<0.7`$, this soliton family is subject to oscillatory instability which is accounted for by the so-called Hamiltonian Hopf bifurcation . In the latter, the collision of two pairs of eigenvalues of opposite Krein signature  results in the generation of a quartet of genuinely complex eigenvalues. The strongest instability, with $`\mathrm{Im}\omega 0.107`$ occurs at $`V_00.5`$.

Similarly, we examined the evolution of the solutions, setting $`V_0=1`$ and $`\mu =1.5`$ in Eq. (3) and performing continuation in $`ϵ`$, starting from $`ϵ=0`$. In this case too, oscillatory instability was found in a finite interval, $`0.05<ϵ<1.61`$. The maximum instability growth rate (which is very large in this case) is $`0.48`$, and occurs at $`ϵ0.6`$. The results for this case are shown in Fig. 2.

In these cases, the oscillatory instability arises quite naturally. Indeed, similar instabilities had been observed in the discrete counterpart of the system and, in fact, in the next section we will give a theoretical justification for their occurrence. Another interesting numerical finding is that the TLM solutions appear to be much more robust than one would initially suspect. In particular, there seem to exist TLM solutions even when the OL has been almost completely smeared out by the harmonic trap (see e.g., the last panel of Fig. 2). We will return to this point in the following section.

We also performed the continuation of the TLM solutions with respect to the chemical potential $`\mu `$, which has yielded results displayed in Fig. 3. In this case, $`ϵ=0`$ and $`V_0=1`$ were chosen. The solutions exist in the region $`\mu <\mu _{\mathrm{max}}0.125`$. The fact that the cutoff value $`\mu _{\mathrm{max}}`$ is positive is explained by a contribution of the mean value of the OL potential $`V_0\mathrm{cos}^2x`$ in Eq. (3). In this case, oscillatory instability is found in two finite intervals, $`0.375<\mu <0.1`$, and $`0.075<\mu <0`$.

### B Bound states of two TLM solitons

Complexes which may be considered as a bound state of two TLMs, or simply as a concatenation of three fundamental solitons with phase difference $`\pi `$ between adjacent ones, were also studied. In this case, for symmetry purposes, $`\mathrm{cos}x`$ in Eqs. (1) and (3) was replaced by $`\mathrm{sin}x`$ (recall we have set $`\lambda =2\pi `$), and three pulses in the initial ansatz (4) were set at $`\xi _j=\pi ,0,\pi `$. We performed the continuation in $`ϵ`$, which has demonstrated that the bound state persists to very large values of $`ϵ`$, see Fig. 4. The branch becomes unstable, through a quartet of complex eigenfrequencies, at $`ϵ>0.015`$. A second instability, accounted for by another eigenvalue quartet, sets in at $`ϵ>0.075`$. This second oscillatory instability can be easily explained: as we argue below, for each pair of anti-phase pulses there exists a pair of eigenvalues with negative Krein signature, which is known to give rise to a Hamiltonian Hopf bifurcation . The strongest instability is found at $`ϵ=0.23`$, when the two instability growth rates are $`0.39`$ and $`0.35`$. The first quartet returns to stability at $`ϵ>0.29`$, while the second quartet is stabilized for $`ϵ>0.33`$. The bound states are completely stable for $`ϵ>0.33`$.

### C Evolution of unstable TLM solitons

In the cases for which the TLM solitons or their bound states were found to be unstable (due to the Hamiltonian Hopf bifurcations), we simulated nonlinear development of the instability in the framework of the full GP equation (1). To this end, a finite-difference scheme with $`dx=0.2`$ \[the same stepsize as in dealing with the stationary equation (3)\] was used, and time integration was performed by means of the fourth-order Runge-Kutta scheme with a time step $`dt=0.001`$. Three examples of the evolution of unstable configurations are shown in Fig. 5. The two top ones correspond, respectively, to the strong and weak oscillatory instability of the TLM soliton. In both of these cases, the TLM state evolves into a single fundamental soliton. A difference between the two cases is that, in the former one, with the relatively weak OL, the newly formed soliton keeps moving between two minima of the potential, while in the latter case , which corresponds to a stronger OL, the fundamental soliton becomes pinned. Finally, the bottom subplot of Fig. 5 demonstrates that the oscillatory instability of a bound state of two TLM solitons transforms it into a single fundamental soliton with a large amplitude.

## III Analytical results

### A The minimum size of the TLM state

The existence of the TLM state in the GP equation with OL can be explained, by means of perturbation theory, as a bound state of two fundamental solitons with phase difference $`\pi `$ between them. Here, we consider the case of $`ϵ=0`$ (no parabolic trap), but the analysis can be easily generalized to include the trap. If the present model is considered as a perturbed NLS equation, we can use an ansatz in the form of Eq. (4) for two solitons. The system’s Hamiltonian reads

$`H={\displaystyle _{\mathrm{}}^{\mathrm{}}}\left[|u_x|^2{\displaystyle \frac{1}{2}}|u|^4+V_0\mathrm{cos}^2\left({\displaystyle \frac{2\pi x}{\lambda }}\right)|u|^2\right]𝑑x.`$

Recall that we have set $`ϵ=0`$; here, for generality’s sake, we do not fix $`\lambda =2\pi `$. Then, known methods (see a review in Ref. ) yield an effective net potential of the interaction between two solitons and interaction of both solitons with the OL. This potential is realized as a part of the Hamiltonian that depends on the coordinates $`\xi _1`$ and $`\xi _2`$ of the solitons and phase difference $`\mathrm{\Delta }\varphi `$ between them (cf. a similar potential for solitons in the discrete-NLS model found in Ref. ) and is of the form:

$`V_{\mathrm{eff}}(\xi _1,\xi _2;\mathrm{\Delta }\varphi )`$ $`=`$ $`8\eta ^3\mathrm{cos}\left(\mathrm{\Delta }\varphi \right)\mathrm{exp}\left({\displaystyle \frac{\eta }{\sqrt{2}}}\left|\xi _1\xi _2\right|\right)`$ (7)
$`+{\displaystyle \frac{2\sqrt{2}\pi ^2V_0}{\lambda \mathrm{sinh}\left(\frac{2\sqrt{2}\pi ^2}{\lambda \eta }\right)}}\left[\mathrm{cos}\left({\displaystyle \frac{4\sqrt{2}\pi }{\lambda }}\xi _1\right)+\mathrm{cos}\left({\displaystyle \frac{4\sqrt{2}\pi }{\lambda }}\xi _2\right)\right].`$
The second term in this expression is the Peierls-Nabarro potential (see, e.g., Ref.  and references therein) induced by the OL. Equilibrium positions of the two-soliton system are determined by equations

$$\frac{V_{\mathrm{eff}}}{\xi _1}=\frac{V_{\mathrm{eff}}}{\xi _2}=\frac{V_{\mathrm{eff}}}{\left(\mathrm{\Delta }\varphi \right)}=0.$$
(8)
Straightforward consideration of Eqs. (8) and (7) shows that a bound state of two solitons with $`\mathrm{\Delta }\varphi =\pi `$, which corresponds to a TLM-like pattern, may exist if the separation $`L\left|\xi _1\xi _2\right|`$ between them exceeds the minimum value,

$$L_{\mathrm{min}}=\frac{\sqrt{2}}{\eta }\mathrm{ln}\left[\frac{\lambda ^2\eta ^4}{\sqrt{2}\pi ^3V_0}\mathrm{sinh}\left(\frac{2\sqrt{2}\pi ^2}{\lambda \eta }\right)\right].$$
(9)
This prediction for the minimum separation between fundamental-soliton components of the TLM state was tested numerically. As a typical example, in Fig. 6 we show the variation of $`L_{\mathrm{min}}`$ with $`V_0`$ for $`\lambda =2\pi `$ and $`\eta =\sqrt{3}`$. The distance between the centers of the solitons was extracted from numerical data as the distance between two local maxima of $`\left|u(x)\right|`$, an inherent discretization error of $`\pm 0.2`$ being imposed by the finite-difference scheme. As is seen, the agreement between the theoretical prediction and the numerically approximated values of $`L_{\mathrm{min}}`$ is quite good.

A more detailed analytical consideration of dynamical properties of a perturbed bound state (following the lines of Ref. ) shows that, within the framework of the present approximation, the bound state is stable. Moreover, the existence of stable bound states of three fundamental solitons, which correspond to the bound state of two TLMs considered above, can also be demonstrated by means of this approach.

### B Oscillatory instabilities of the single-TLM and double-TLM states

It is well known (see e.g., Ref.  and references therein) that the linearized equations for the perturbations defined in Eq. (5) can be rewritten, in terms of new variables $`U=a+b^{}`$ and $`W=ab^{}`$, as

$`\omega U`$ $`=`$ $`L_{}WW_{xx}v^2W+(V(x)\mu )W,`$ (10)
$`\omega W`$ $`=`$ $`L_+UU_{xx}3v^2U+(V(x)\mu )U_n.`$ (11)
The Krein signature of the eigenvalue can be written as

$`K=\mathrm{sign}\left({\displaystyle _{\mathrm{}}^{\mathrm{}}}(UW^{}+WU^{})𝑑x\right).`$

Analytical consideration is possible for $`\omega `$’s that are real; then $`U`$ and $`W`$ are also real, hence $`K`$ amounts to the expression

$$K=\mathrm{sign}\left(_{\mathrm{}}^{\mathrm{}}UW𝑑x\right).$$
(12)
As was done above, we approximate a TLM solution as a superposition, $`v(x)=P(x)+N(x)`$, of two stationary fundamental-soliton waveforms, positive $`P(x)`$ and negative $`N(x)`$, the underlying assumption being a weak overlap between $`P(x)`$ and $`N(x)`$. Since $`P`$ and $`N`$ are, essentially, individual solitons, in the lowest-order approximation they obey equations $`\left[V(x)\mu \right]P=P_{xx}+P^3`$ and $`\left[V(x)\mu \right]N=N_{xx}+N^3`$. Summing them up, and taking into account that the product $`\left|PN\right|`$ is everywhere much smaller than $`P^2+N^2`$ (due to the assumed weak overlap between $`P`$ and $`N`$), we conclude that, up to the same lowest-order accuracy, the function $`P(x)+N(x)`$ satisfies the equation $`L_{}(P+N)=0`$ \[recall that the operator $`L_{}`$ is defined in Eq. (10)\]. Thus, in the lowest approximation, $`P(x)+N(x)`$ is a zero mode of the operator $`L_{}`$. In a similar way, the following approximate result can be obtained for the function $`P(x)N(x)`$:

$$L_{}(PN)=PN(PN).$$
(13)
The linearization around each of the TLM-constituent pulses (if considered in isolation) carries a pair of zero eigenvalues due to the phase invariance. The corresponding eigenfunction, $`W`$, has the shape of the pulse itself ($`P`$ and $`N`$, respectively). There is also a generalized eigenfunction  for each pair. We expect to have two pairs of TLM eigenvalues originating from these pairs of the individual-soliton’s ones, the corresponding eigenfunctions being (in the lowest approximation) $`W=P+N`$ and $`W=PN`$. The former one is related to the phase invariance of the system (in other words, to the conservation of the total number of atoms), which means that the corresponding eigenvalue pair remains exactly equal to zero, while the other one may become different than zero. Then, it follows from Eq. (13) and Eq. (10) that, in the lowest approximation, $`U=PN(PN)`$. In turn, Eq. (12) yields

$$K=\mathrm{sign}\left[_{\mathrm{}}^{\mathrm{}}PN(PN)^2𝑑x\right].$$
(14)
Since, by definition, the signs of $`P(x)`$ and $`N(x)`$ are opposite, Eq. (14) implies that the Krein signature of this phase eigenmode is negative. A known consequence of this  is that an oscillatory instability will be generated by collision of this nonzero eigenvalue with other discrete ones (and/or those belonging to the continuous spectrum) that carry a positive Krein signature.

The above consideration indicates that, for every pair of the interacting $`\pi `$-out-of-phase pulses, there exists an eigenvalue with negative Krein signature and, consequently, a potentially ensuing Hamiltonian Hopf bifurcation. This conclusion fully agrees with our numerical results presented above, that have identified a single oscillatory instability for the single TLM soliton, and two such instabilities in the case of bound states of two TLM solitons.

### C TLM solitons in the absence of the OL

Finally, it is possible to present analytical arguments justifying the existence of TLM solitons (both single and multiple ones) even when the harmonic trap is overwhelmingly stronger than the OL potential (see the numerical results shown in the last panels in Figs. 2 and 4), or when the OL potential is simply absent. Within the framework of the same perturbation theory which was employed above to derive Eq. (7), each of the pulses constituting the TLM state is then subject to the action of two forces. One of them is exerted by the magnetic trap, while the other comes from the interaction between the pulses. Accordingly, the effective potential is

$`V_{\mathrm{eff}}(\xi _1,\xi _2;\mathrm{\Delta }\varphi )`$ $`=`$ $`8\eta ^3\mathrm{cos}\left(\mathrm{\Delta }\varphi \right)\mathrm{exp}\left({\displaystyle \frac{\eta }{\sqrt{2}}}\left|\xi _1\xi _2\right|\right)`$ (16)
$`+4\sqrt{3}ϵ\left(\xi _1^2+\xi _2^2\right).`$
Predictions following from this approximation were tested against direct numerical results obtained in the absence of the OL, i.e., with $`V_0=0`$ in Eq. (1). In Fig. 7, an example is displayed for $`ϵ=0.2`$ and $`\mu =1.5`$. The top panel shows the effective potential (16), which, in this case, has a minimum at $`\xi _{1,2}=\pm 1.35`$. The middle and bottom panels present a numerical solution, with maxima of $`\left|u(x)\right|`$ at $`\xi _{1,2}\pm 1.4`$. Thus, the perturbative approximation not only explains the existence of the TLM state in attractive BECs even in the absence of OL, but also accurately predicts the location of centers of the pulses whose concatenation generates the TLM structure. This, in turn, supports the numerical findings presented in section II, which suggested that TLM solitons in the attractive BEC would indeed persist even for very tight magnetic traps. Thus, we conclude that the TLM solitons are a very robust feature of the system.

## IV Conclusions

In this work we have demonstrated that counterparts of twisted localized modes (TLMs), which were earlier found in the discrete one-dimensional NLS equation, can exist as robust objects in attractive, effectively one-dimensional Bose-Einstein condensates (BECs), provided that the condensate is confined by the optical lattice (OL) and/or parabolic magnetic trap. The TLMs in BECs may be considered as bound states of two fundamental solitons with a phase difference $`\pi `$ between them, or, alternatively, as a dark soliton embedded in a broader bright one. The family of the TLM solitons and their stability were investigated in detail, by means of numerical continuation in the strengths of the OL or magnetic-trap parameters, as well as in the chemical potential. Stable bound states of two TLMs, alias a bound state of three fundamental solitons with a $`\pi `$ phase shift between adjacent ones, were also found. In the cases when the TLMs are unstable, the development of the instability was monitored by means of direct simulations, with a conclusion that the TLM rearranges itself into a fundamental soliton.

We have also developed an analytical approach, which treats each fundamental soliton as a quasi-particle subject to the action of effective forces exerted by the OL, magnetic trap, and interaction with another soliton. The analysis predicts the existence of the stable bound state of two fundamental solitons with a phase shift of $`\pi `$ between them. The smallest possible distance between the bound solitons, predicted by this approach, is in good agreement with the numerically found size of the corresponding TLM. The analysis can be readily extended to predict that multi-pulse bound states with the phase shift $`\pi `$ between adjacent solitons also exist and may be stable.

The origin of the instability of the TLMs, in the case when they are unstable, was also examined by means of an analytical approximation, and was found to stem from Hamiltonian Hopf bifurcations, which occur when eigenvalues with negative Krein signature collide with positive-signature ones. This result was based on the existence of a negative-Krein-sign eigenvalue set, that was constructed by means of the perturbative analysis.

It will be very interesting to generalize these considerations to higher dimensions – in particular to two-dimensional settings. In that case, a relevant issue is the possible existence of a 2D counterpart of the TLM in the form of a bright vortex (which was found in the 2D discrete NLS equation, where it was shown to be stable under certain conditions ). A very recent preliminary result is that bright vortices, stabilized by a 2D optical lattice, can indeed be found .

## Acknowledgements

One of the authors (B.A.M.) appreciates valuable discussions with B. Baizakov and M. Salerno, and hospitality of the Department of Physics at Universitá di Salerno (Italy). The work of this author was supported, in a part by a Window-on-Science grant from the European Office of Aerospace Research and Development (EOARD) of US Air Force. This project was also supported by a University of Massachusetts Faculty Research Grant and NSF-DMS-0204585 (P.G.K.), the Special Research Account of the University of Athens (D.J.F.), the Binational (US-Israel) Science Foundation, under grant No. 1999459 (B.A.M.) and by the AFOSR Dynamics and Control (I.G.K.). Work at Los Alamos is supported by the US DoE.


text:
# Erratum: Superconductivity in CoO2 Layers and the Resonating Valence Bond Mean Field Theory of the Triangular Lattice t-J Model [Phys. Rev. B 68 , 104508 (2003)]

Our recent paper  contains an error in the formula for the high frequency Hall constant for the experimentally relevant case of electron doping. Eq(3) giving the high frequency Hall constant for the case of electron filling, i.e. electron density $`n>1`$, should correctly read

$$R_H^{}=\frac{v}{4|e|}\frac{k_BT}{t}\frac{1+\delta }{\delta (1\delta )}\text{(Cases A(i) or B(ii))}\text{Eq(3’)}$$

where $`\delta =|1n|`$. This differs essentially from the quoted answer by a minus sign that was missed by us, and by a factor of 2. This equation is obtained from Eq(2) in the paper, calculated for the case of hole doping $`n<1`$, by a particle hole transformation. Performing a canonical transformation $`c_{i,\sigma }c_{i,\sigma }^{}`$ in the Kubo formulas for models of lattice fermions (see e.g.  Eqn(1,2) ) yields:

$$R_H(n,t,\omega ,T)=R_H(2n,t,\omega ,T).$$

The crucial overall minus sign arises from the fact that the Hall constant involves a three current propagator, and each current picks up a minus sign under the transformation. Incidentally this identity results in a vanishing at half filling for bipartite lattices of the Hall function $`R_H(1,t,\omega ,T)`$ for any value of $`T,\omega `$ and the interaction strength.

The changed sign in Eq(3’) has a few implications which we would like to state explicitly. The recent experiment of Wang, Rogado, Cava and Ong on the transport Hall constant, actually sees the predicted linear T behaviour over a wide temperature range $`200^0KT400^0K`$ with a positive slope. From Eq(3’) above, we deduce that the appropriate $`tJ`$ model description of the system has $`t<0`$, i.e. belongs to case B(ii) in our notation, at least at the composition $`x.7`$.

We note that with this correction, the sign of the hopping $`t`$ suggested by our calculation is the same as the one favoured by Refs, at least for the non superconducting range $`x.7`$. In this composition range, photoemission data  also suggest the same sign of $`t`$, with an unoccupied region around the center of Brillouin zone, namely the $`\mathrm{\Gamma }`$ point.


text:
# Electron spin precession in semiconductor quantum wires with Rashba spin-orbit coupling

## 1 Introduction

The control of the Rashba spin-orbit coupling is one of the most promising tools for the manipulation of the spin of the carriers within semiconductor heterostructures. In the last years, there have appeared several proposals of spin-based devices relying on this mechanism peet02 ; egue02 . Among them we cite the spin-FET, first proposed by Datta and Das datt90 , the spin filters koga02  and the spin guides lolo02 .

The origin of the Rashba spin-orbit coupling in III-V semiconductor heterostructures lies in their asymmetry. The lack of space symmetry causes a local electric field perpendicular to the plane of the heterostructure. The associated relativistic correction makes the electrons feel an effective magnetic field which separates in energy the different spin states ras60 . This dependence of the Rashba spin-orbit coupling on the interface electric field remains as a great advantage for spin control in semiconductors. Indeed, the intensity of the coupling can be controlled by applying a vertical electric field to the heterostructure, as proved experimentally by Nitta et al. nitt97 .

In this work, we shall focus on the spin precessional properties of conduction electrons in a quantum wire having a finite width. The problem of the quantum wire with Rashba coupling has already been studied by Moroz and Barnes moro99 , who addressed the ballistic conductance and the subband structure of the wire; and by Mireles and Kirczenow mire01 , who analyzed the spin-dependent transport within a tight binding model. Following a different approach, Governale and Zülicke gove02  have also treated the problem, and, more recently, Egues et al egue02b  have proposed to take advantadge of the coupling of the first two transverse subbbands to create a spin-FET with enhanced spin control. In this work, we shall concentrate on how the coupling of the different transverse subbands, due to the Rashba term, affects the electron spin precession.

The Hamiltonian representing the Rashba spin-orbit coupling in terms of the electron momentum ($`𝐩`$) and the Pauli matrices ($`\sigma `$’s) reads

$$_{SO}=\frac{\lambda _R}{\mathrm{}}\left(p_y\sigma _xp_x\sigma _y\right),$$
(1)
where $`\lambda _R`$ represents the intensity of the spin-orbit coupling, depending on the heterostructure’s vertical electric field. To stress the experimental feasibility we have assumed a value $`\lambda _R=1.03\times 10^9`$ eV cm, as reported in Ref. nitt97  for an In0.53Ga0.47As/In0.52Al0.48As quantum well. Other experimental parameters corresponding to this sample are: $`m^{}=0.05m_e`$, $`ϵ=13.9`$ for InGaAs, and the Fermi wavevector $`k_F=3.5\times 10^6`$ cm-1, corresponding to a 2DEG density $`n_s2\times 10^{12}`$ cm-2. All calculations discussed below have been obtained using this parameter set.

We assume a wire oriented along the $`y`$ axis with confinement achieved in the $`x`$ direction. In the simplified model of a one-dimensional wire the terms containing $`x`$ or $`p_x`$ are dropped and the remaining Hamiltonian

$$_{1D}=\frac{p_y^2}{2m^{}}+\frac{\lambda _R}{\mathrm{}}p_y\sigma _x$$
(2)
is compatible with the observables $`\{p_y,\sigma _x\}`$. That is, the energy eigenstates are characterized by their $`x`$ spin component and $`y`$ momentum ($`\mathrm{}k_y`$). In this system, the state of an electron having well-defined momentum and with arbitrary spin orientation decomposes into a combination of the spin eigenstates $`|s_x+`$ and $`|s_x`$, leading to a spin precessional motion about the $`x`$ axis, i.e., in the $`S_y`$-$`S_z`$ plane. The spin evolution is characterized by an angular frequency $`\omega _P`$, corresponding to the energy difference between the two spin eigenstates, $`\mathrm{}\omega _P=2\lambda _Rk_y`$. This precession frequency leads to the well-known formula for the angular modulation of spin orientation $`\theta _R=2m^{}\lambda _RL/\mathrm{}^2`$, for a channel length $`L`$.

When a two-dimensional free motion of the electrons is considered, the eigenstates of the Hamiltonian

$$_{2D}=\frac{\left(p_x^2+p_y^2\right)}{2m^{}}+\frac{\lambda _R}{\mathrm{}}\left(p_y\sigma _xp_x\sigma _y\right)$$
(3)
have again well-defined momentum $`\mathrm{}𝐤`$ and spin orientation, but the latter depends now on the wavevector $`𝐤`$. Using spinorial notation in the usual basis of $`|+`$ and $`|`$ states for the $`S_z`$ spin component, the wave functions and energies for a given $`𝐤`$ read

$`\chi _{𝐤\pm }`$ $`=`$ $`e^{i𝐤𝐫}\left(\begin{array}{c}1\\ ie^{i\varphi }\end{array}\right),`$ (6)
$`\epsilon _{𝐤\pm }`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{2m^{}}}(k_x^2+k_y^2)\pm \lambda _R|𝐤|,`$ (7)
where $`\varphi =\mathrm{atan}(k_y/k_x)`$ is the polar angle for $`𝐤`$. Note that the spinors (6) are eigenstates of a spin component orthogonal to $`𝐤`$, i.e., the spin orientation depends on the spatial wavevector, always pointing in perpendicular direction to $`𝐤`$. As in the 1D case, a state having well-defined momentum $`\mathrm{}𝐤`$ and arbitrary spin orientation decomposes into a linear combination of the above 2D eigenspinors $`\chi _{𝐤\pm }`$. Therefore, the evolution under $`_{2D}`$ leads again to a spin precession of frequency $`\omega _P`$ given by the energy difference between the corresponding ’up’ and ’down’ eigenstates. Namely, $`\mathrm{}\omega _P=2\lambda _R|𝐤|`$; the same result of the 1D case. An important difference is, however, that the plane of spin precession depends in 2D on the spatial orientation of $`𝐤`$. Having analyzed the bulk limits, the question that arises is how the above pure spin precession is modified when the complete translational invariance of the homogeneous 1D and 2D systems is lost in a wire with a finite width.

## 2 The Model

In order to model the quantum wire we maintain the translational invariance in the longitudinal coordinate ($`y`$) while, in the transversal direction ($`x`$) this symmetry is broken by a parabolic confinement

$$V(x)=\frac{1}{2}m^{}\omega _0^2x^2,$$
(8)
where $`\omega _0`$ determines the width of the wire. The Hamiltonian representing the quantum wire reads

$$=\frac{𝐩^2}{2m^{}}+\frac{1}{2}m^{}\omega _0^2x^2+\frac{\lambda _R}{\mathrm{}}\left(p_y\sigma _xp_x\sigma _y\right).$$
(9)
As our system still retains one spatial symmetry we can reduce the spinorial 2D problem to a 1D one by factorizing out the $`y`$ dependence,

$$\left(\begin{array}{c}\psi _{nk_y}(𝐫)\\ \psi _{nk_y}(𝐫)\end{array}\right)=e^{ik_yy}\left(\begin{array}{c}\varphi _{nk_y}(x)\\ \varphi _{nk_y}(x)\end{array}\right).$$
(10)
The continuum index $`k_y`$ in Eq. (10) labels the state of free motion in longitudinal direction while the discrete index $`n`$ accounts for the different transversal subbands of the confining dimension. Introducing this general form of energy eigenspinor we obtain the reduced Hamiltonian for the spinorial transverse dependence.

$`_{\mathrm{𝑡𝑟}}`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2k_y^2}{2m^{}}}+{\displaystyle \frac{p_x^2}{2m^{}}}+{\displaystyle \frac{m^{}\omega _0^2}{2}}x^2`$ (11)
$`+`$ $`{\displaystyle \frac{\lambda _R}{\mathrm{}}}\left(\mathrm{}k_y\sigma _xp_x\sigma _y\right)`$
The transverse Hamiltonian, Eq. (11), can be analyzed by seeking direct solutions to the stationary Schrödinger equation or, equivalently, solving the time-dependent Schrödinger equation for given initial conditions. We shall use the latter approach since it yields a suitable scenario to address the dynamical spin evolution. For a given wavevector $`k_y`$ the coupled time-dependent equations for the two spin components of any spinorial wavepacket: $`f_{k_y}(x)`$ and $`f_{k_y}(x)`$ read

$`i\mathrm{}_t\left(\begin{array}{c}f_{k_y}(x,t)\\ f_{k_y}(x,t)\end{array}\right)=`$ (14)
$`\left(\begin{array}{cc}_0^{(x)}& \lambda _R\left(k_y+_x\right)\\ \lambda _R\left(k_y_x\right)& _0^{(x)}\end{array}\right)\left(\begin{array}{c}f_{k_y}(x,t)\\ f_{k_y}(x,t)\end{array}\right),`$ (19)
where

$$_0^{(x)}=\frac{\mathrm{}^2}{2m^{}}_x^2+\frac{m^{}\omega _0^2}{2}x^2.$$
(20)
Note that in $`_0^{(x)}`$ we have dropped the longitudinal kinetic contribution, not important for time evolution.

It is well known that in quantum mechanics time evolution can provide complete information about the stationary states of a Hamiltonian through spectral analysis. Indeed, for a time-independent Hamiltonian an arbitrary wavepacket can be expanded as a series in the stationary energy eigenstates

$$\left(\begin{array}{c}f_{k_y}(x,t)\\ f_{k_y}(x,t)\end{array}\right)=\underset{n}{}A_{nk_y}e^{\frac{i}{\mathrm{}}\epsilon _{nk_y}t}\left(\begin{array}{c}\varphi _{nk_y}(x)\\ \varphi _{nk_y}(x)\end{array}\right).$$
(21)
Using the fast Fourier transform (FFT) we can extract the harmonic frequencies ($`\epsilon _{nk_y}`$ energies) from the evolution of an arbitrary wavepacket. Furthermore, if this analysis is repeated at different spatial points the space components of the energy eigenspinor $`\varphi _{nk_y}(x)`$ and $`\varphi _{nk_y}(x)`$, can also be obtained.

## 3 Results

The Hamiltonian of Eq. (9), apart from kinetic and potential terms, contains the longitudinal and transversal spin-orbit contributions coupling $`x`$ and $`y`$ spatial motions with the electron spin. It is worth noticing that there is an analytical limit to the solution of Eq. (19) when $`p_xf_{k_y}\mathrm{}k_yf_{k_y}`$. In this case the solutions are similar to those discussed above for $`_{1D}`$, i.e., the eigenstates are still plane waves in $`y`$ direction with spin oriented along $`x`$ but now with an harmonic oscillator transverse profile. The corresponding energy spectrum is also analytically known

$$\epsilon _{nk_ys}^{(0)}=\frac{\mathrm{}^2k_y^2}{2m^{}}+\left(n+\frac{1}{2}\right)\mathrm{}\omega _0+\lambda _Rk_ys,$$
(22)
where $`s=+1`$ for $`+x`$ spin orientation (up) and $`s=1`$ for $`x`$ (down). If we substract the longitudinal kinetic term $`\mathrm{}^2k_y^2/2m^{}`$, common to both spin orientations, the remaining band structure consists of doubly split transverse subbands, linear in $`k_y`$ and having a common origin at $`(n+1/2)\mathrm{}\omega _0`$. Another analytical limit can be obtained neglecting $`\lambda _Rp_y\sigma _x`$ in (9); in this case, spin-orbit interaction only produces a little constant shift in the energy levels that does not depend on spin orientation.

As a check of our numerical procedure, Fig. 1 shows the results for a simulation using the above simplified Hamiltonian (neglecting $`\lambda _Rp_x\sigma _y`$ in Eq. (9)). The parameters used are: $`\mathrm{}\omega _0=3.5`$ meV, $`k_y=1.0\times 10^6`$ cm-1 and the above mentioned spin-orbit intensity. The initial input for the time-evolution is a Gaussian-shaped wavepacket having spin oriented in $`+z`$,

$$\left(\begin{array}{c}f_{k_y}(x,0)\\ f_{k_y}(x,0)\end{array}\right)\frac{1}{2\sigma \sqrt{\pi }}e^{\frac{(xx_0)^2}{2\sigma ^2}}\left(\begin{array}{c}1\\ 0\end{array}\right).$$
(23)
Using a finite $`x_0`$ and $`\sigma \sqrt{\frac{\mathrm{}}{m^{}\omega _0}}`$ we ensure that the initial wavepacket is composed of several transversal eigenstates in the low energy region. Figure 1 shows the spectrum corresponding to the wavepacket’s oscillations at an arbitrary $`x`$ point, as well as the transversal densities extracted for the two lowest peaks. An excellent agreement between numerical and analytical energies and eigenspinors is found. Figure 2 represents the dynamical spin evolution for the above initial spinor, corresponding to a pure precession in excellent agreement with the analytical frequency $`\mathrm{}\omega _P=2\lambda _Rk_y`$.

Note that, in this analytical regime, spin precession is independent of the transversal state, as the energy gaps are common to the different subbands. This statement implies that spin evolution is independent of the transversal profile of the injected particles. Numerically, we have checked that arbitrary initial profiles for the $`x`$ dependence lead to the same spin dynamics.

The above analytical regime no longer holds when the transverse motion couples with the spin through the term $`\lambda _Rp_x\sigma _y`$ in Eq. (9). Complete analytical solutions for the Hamiltonian including the full spin-orbit interaction are not available and numerical calculations are needed to solve the problem. The first question arising is how the band structure of the simplified Hamiltonian is modified when transversal spin-orbit coupling is included. Figure 3 shows the corresponding band structure for the same wire parameters considered above. In order to clarify the effect of the transversal coupling the longitudinal kinetic term $`\mathrm{}^2k_y^2/2m^{}`$ has been subtracted and the analytical energies are shown with dotted lines. It can be seen that for low values of $`k_y`$ there is good agreement between the simplified and full Hamiltonians; but this is not true when we approach a region of band crossings. When the $`(n+1)`$-th down subband reaches the up branch of the $`n`$-th subband, the transversal coupling produces a large effect, removing the band degeneracy and opening a gap between these two branches, $`(n+1,)`$ and $`(n,)`$ moro99 . The gap opened clearly depends on the transversal subbands involved; the higher transversal state the larger the gap, as easily noticed from Fig. 3. It can also be observed that the range of $`k_y`$’s where there is a good agreement between full calculations and analytical model disminishes as the subband index increases. These changes on the level structure suggest significant modifications of the spin precession due to the coupling with transversal motion.

As shown by Fig. 3, transversal coupling effects maximize at the analytical model crossings, occurring approximately when $`2\lambda _Rk_y=\mathrm{}\omega _0`$. The fulfillment of this condition establishes a rough criterion to quantify the relevance of the transversal spin-orbit coupling. Nevertheless, it must be also stressed that the coupling depends sensitively on the subband index $`n`$.

The influence of the transversal coupling can also be observed in the structure of the eigenspinors. In the analytical model the eigenfunctions are states having well-defined spin orientation. The transversal spin-orbit coupling mixes the different subbands leading to a more complex structure of the eigenstates. Figure 4 represents the densities corresponding to three lowest eigenstates in the region of maximum band anticrossing, i.e., $`k_y=1.70\times 10^6`$ cm-1. The lowest energy state corresponds to a branch without crossings and, therefore, its eigenvalue and eigenspinor are very similar to those of the analytical model; spin density mainly oriented in $`x`$ direction with a residual structure in $`z`$ orientation. The other two eigenspinors differ substantially from those of the analytical model; spin has no predominant orientation and displays a complex texture very similar to those found in Ref. gove02 . Charge and spin densities are related to the eigenspinor components through the usual expressions

$`\rho (x)`$ $`=`$ $`\varphi _{nk_y}(x)^2+\varphi _{nk_y}(x)^2`$
$`m_z(x)`$ $`=`$ $`\varphi _{nk_y}(x)^2\varphi _{nk_y}(x)^2`$
$`m_x(x)`$ $`=`$ $`2\mathrm{R}\mathrm{e}\left\{\varphi _{nk_y}(x)\varphi _{nk_y}^{}(x)\right\}.`$ (24)
Until now we have presented the effects produced by the full spin-orbit coupling on the level scheme with respect to the analytical model, having results in good agreement with those of Ref. moro99 . The question that we address next is how this modified band structure reflects on the spin dynamics. We consider first the case of a particle having well-defined spin orientation at initial time and $`y`$-momentum in the range of band anticrossings, i.e., an spinor as in (21) evolving with $`k_y=1.70\times 10^6`$ cm-1 in Eq. (19). The initial condition thus reads

$$f_{k_y}(x,0)\frac{1}{\pi ^{1/4}}\left[\frac{m^{}\omega _0}{\mathrm{}}\right]^{1/4}e^{\frac{m^{}\omega _0}{2\mathrm{}}x^2}\left(\begin{array}{c}1\\ 0\end{array}\right).$$
(25)
The spatial part of this injected spinor is equal to that of the lowest analytic subband. In this way, we ensure that it is in the low energy range, very close to the first transversal subband.

The dynamical evolution of this injected particle is shown in Fig. 5. The spin trajectory on the $`S_y`$-$`S_z`$ plane is quite involved, manifesting a multifrequential evolution. Indeed, the Fourier transform of the series shows that spin evolves in time as a superposition of several frequencies corresponding to the energy transitions between the different anticrossed branches. Since the gap opened by the transversal spin-orbit coupling varies from one subband to another, several precession frequencies are found in the anticrossing region. On the contrary, the analytical model predicts a common transition energy for all the transversal subbands (see Fig. 3), leading to a unique precession frequency. This feature proves that spin evolution in 2D wires noticeably depends on the transversal state of the particle. In the particular case of Fig. 5 the precession is characterized by two dominant frequencies, corresponding to the transitions between the lowest branch and the next two higher branches, represented by dark-gray arrows in Fig. 3. This two-frequency evolution yields a variable spin amplitude, in agreement with the displayed trajectory. The spin-orbit coupling makes the precessional frequencies depend on $`k_y`$ in a non-trivial way. At the same time, it gives rise to a $`k_y`$-dependent spin modulation angle that would lead to a certain degree of decoherence in the spin transport through the wire.

Figure 6 shows the precessional spectra corresponding to an initial spinor having well-defined spin (in $`+z`$ direction) and different transversal states. The precessional spectrum corresponding to the first transversal subband ($`n=0`$) is that represented in Fig. 5. The precessional spectra for the higher transversal states, $`n=1`$ and $`n=2`$, exhibit an enhanced multi-frequential character, since more transitions significantly contribute to the spin evolution. It can also be seen that the strength spreads over a higher frequency interval as the subband index $`n`$ is increased. For low $`k_y`$’s, out of the anticrossing region, the spin evolution recovers a single-frequency behaviour, with slight deviations from the analytical results (cf. Fig. 2). Note that, in general, a pure up ($`+z`$) spin injected into the wire will couple with the different spinor eigenstates and, therefore, the multifrequential evolution can not be avoided.

Another question of relevance is how the transversal size influences the precession. From Fig. 4 we can estimate the width of the wire in $`90120`$ nm for $`\mathrm{}\omega _0=3.5`$ meV. If we reduce the value of $`\omega _0`$ and maintain the intensity of spin-orbit coupling the crossings of transversal bands become more frequent as subbands are less spaced. Figure 7 shows the corresponding band structure in the case $`\mathrm{}\omega _0=0.70`$ meV; from the densities of the eigenspinors an estimated width of $`200250`$ nm can be obtained for this wire. Deviations from the analytical band structure are clearly visible in spite of the reduced $`k_y`$ range displayed. This feature is reflected in the precessional spectra which are fragmented even for relatively low $`k_y`$’s. By contrast, a wire of $`\mathrm{}\omega _0=14`$ meV, having an estimated width of $`4060`$ nm, shows a band structure very close to that of the analytical decoupled Hamiltonian, as displayed in Fig. 8. This also reflects on the precessional spectra, which are single-moded and match the value given by the analytical formula $`\mathrm{}\omega _P=2\lambda _Rk_y`$.

In practical implementations of the Datta and Das spin transistor datt90  the above discussed spin decoherence at subband crossings would strongly affect the device operation and, therefore, it should be carefully considered in the design of appropriate working parameters. We can quantify the decoherence by representing the spin amplitude after the first spin flip, i.e., after the first $`\pi `$ rotation. Figure 9 represents this quantity as a function of wire width ($`\omega _0`$) for a $`k_y`$ value in the anticrossing region, $`\mathrm{}\omega _0/2\lambda _R`$. It can be seen that as the wire becomes wider (lower $`\omega _0`$), the multi-frequential precession leads to an important reduction of the $`z`$ spin after the first inversion.

## 4 Summary

In summary, we have studied the electronic spin precession in quantum wires with Rashba spin-orbit coupling by means of real time simulations. A detailed study of the band structure for this system has been presented and compared with a decoupled analytical model. It has been shown that deviations from the simple spin precession regime are determined by the relative importance of the longitudinal spin-orbit energy ($`\lambda _Rk_y`$) and the transversal confining energy ($`\mathrm{}\omega _0`$) and that the spin precession fragments into several frequencies whose weights depend on the transversal state of the particle. The dependence of the spin precession on the wire width has also been studied and illustrated through the band structure.

This work was supported by Grant No. BFM2002-03241 from DGI (Spain).