format: text

text:
# 1 Introduction

## 1 Introduction

Pulse propagation in nonlinear birefringent optical fibers is presently an active area of research \[1-12\].Even “single-mode” fibers support two degenerate modes polarized in two orthogonal directions. The degeneracy can be broken introducing deliberately a large amount of birefringence through design modifications giving place to the so called polarization-preserving fibers. Weak light pulses propagating in these media maintain its state of polarization when they are initially polarized along either of the principal axes. However at high intensities nonlinear birefringence can change drastically the propagation dynamics. Birefringence can also be a residual effect due to imperfections in the manufacturing process.

The equations that govern pulse propagation in nonlinear birefringent pulses have been derived by Menyuk . Two main cases of birefringence: i) high and ii) low, have been usually considered separately. In the approximation of high birefringence , phase and group velocities of both field components are considered to be different. At the same time rapidly oscillating nonlinear terms are neglected. On the other hand, the approximation of low birefringence takes into account the difference in phase velocities between the two linearly polarized components, but neglects their difference in group velocities. The full polarization dynamics of solitons in this last approximation has been studied in .

The general case has been analyzed in just a few papers . New vectorial stationary solitons for this case were first discovered in , and later on the whole family of stationary pulselike solutions was numerically determined by Torner et al. . These families of solitons correspond to mutually trapped pulses that propagate with a common group velocity and which therefore can be called walking vector solitons . Stationary solutions as singular points of this infinite dimensional dynamical system play a pivotal role to understand the propagation dynamics of any arbitrary input pulse. Here we generalize the method described in  to find accurate analytical approximations of the solitary wave solutions for the case of high birefringence. We generalize it for the more general case, when no terms are neglected. Our approximate results are compared with existing exact ones  to prove the validity of our approach. We then examine the behavior of the two pulse components when the propagation constant ($`q`$) or the velocity ($`v`$) change. We find that the width of the pulses is determined by $`q`$ whereas the amplitudes ratio $`a_1/a_2`$ is determined by $`v`$ when the energy is higher than a certain value. By means of numerical propagation of the stationary solutions found in this way, we observe that most of them are stable.

The rest of the paper is organized as follows. The problem to be addressed by the variational method is described in section 2. Our variational approach is developed in section 3. In section 4 we present the numerical results and finally in section 5 we briefly summarize the main conclusions.

## 2 Statement of the problem

Pulse propagation in a birefringent optical fiber can be described in terms of two nonlinearly coupled nonlinear Schrödinger equations. In a reference frame traveling along the $`\xi `$ axis with the average group velocity, and in normalized units, this set takes the form

$`iU_\xi =i\delta U_\tau +{\displaystyle \frac{1}{2}}U_{\tau \tau }+\left(|U|^2+A|V|^2\right)U+\beta U+BV^2U^{},`$
$`iV_\xi =i\delta V_\tau +{\displaystyle \frac{1}{2}}V_{\tau \tau }+\left(|V|^2+A|U|^2\right)V\beta V+BU^2V^{},`$ (1)
where $`U`$ and $`V`$ are the slowly varying envelopes of the two linearly components of the field along the $`x`$ and $`y`$ axis respectively, $`\xi `$ is the normalized propagation coordinate, $`\delta `$ is the inverse group velocity difference, $`\beta `$ is half the difference between the propagation constants, $`A`$ is the normalized ratio of the nonlinear susceptibilities, $`B=1A`$, $`\tau `$ is the normalized retarded time, and the asterisk denotes complex conjugation. In linearly birefringent fibers $`A=2/3`$, as we set here. The set of Eqs.2 has at least three integrals of motion , i) the action (total energy), ii) the momentum and iii) the Hamiltonian, which are a consequence of i)the translational invariance of Eqs.(2) relative to phase shifts, ii) invariance in $`\tau `$ and iii) translational invariance in $`\xi `$.

Eqs.2 have two simple linearly polarized pulselike solutions, viz., linearly polarized soliton waves along the slow axis

$$U=\frac{\sqrt{2(q\beta )}}{\mathrm{cosh}[\sqrt{2(q\beta )}(\tau \delta \xi )]}\mathrm{exp}(iq\xi ),V=0$$
(2)
and linearly polarized soliton waves along the fast axis

$$U=0,V=\frac{\sqrt{2(q+\beta )}}{\mathrm{cosh}[\sqrt{2(q+\beta )}(\tau +\delta \xi )]}\mathrm{exp}(iq\xi ).$$
(3)
Their corresponding values for the total energy of the pulse are:

$$Q=_{\mathrm{}}^{\mathrm{}}(|U|^2+|V|^2)𝑑\tau =2\sqrt{2(q\pm \beta )}$$
(4)
The stability of these solutions has been determined in Ref. for $`\delta =0`$, and for arbitrary nonzero values of $`\delta `$ in Ref..

Usually Eqs.2 are written without the last two terms , as terms with coefficient $`\beta `$ can be eliminated by $`z`$ dependent phase transformations of $`U`$ and $`V`$, and this in turn gives the last nonlinear terms rapid phase variations with $`z`$ which in principle allow their neglect. We must however bear in mind that the last terms in Eq.2 are the only ones responsible for energy transfer between both polarizations. And although averaging out the fast oscillatory terms has proven to be a good approximation to describe most of the observed phenomena  in the picosecond regime, we will go further in retaining these terms in the analysis.

## 3 Variational approach

The variational approach when applied to the single Nonlinear Schrödinger Equation (NLS) was introduced by Anderson . Since then it has been widely used for coupled NLS equations. For the specific case of birefringent fibers in the high birefringence approximation it has been used for studying dynamical behaviors , as well as the stationary case .

We look for stationary solutions moving at a common velocity $`v`$ in our frame of reference. They can be written  as

$`U(\tau ^{},\xi )=P_1(\tau ^{})\mathrm{exp}(iq\xi )`$
$`V(\tau ^{},\xi )=P_2(\tau ^{})\mathrm{exp}(iq\xi )`$ (5)
where $`\tau ^{}=\tau v\xi `$ is the common retarded time. Inserting Eqs. 3 into Eqs.2, we get a set of ordinary differential equations (ODE’s) for $`P_{1,2}`$, which reads

$`i(\delta v)\dot{P}_1+(\beta q)P_1+{\displaystyle \frac{1}{2}}\ddot{P}_1+(|P_1|^2+A|P_2|^2)P_1+BP_2^2P_1^{}=0`$
$`i(\delta +v)\dot{P}_2(\beta +q)P_2+{\displaystyle \frac{1}{2}}\ddot{P}_2+(|P_2|^2+A|P_1|^2)P_2+BP_1^2P_2^{}=0`$ (6)
where the overdots indicate derivative respect to $`\tau ^{}`$. Eqs. (3) can be derived (via the Euler-Lagrange equations) from the following Lagrangian $`L=_{\mathrm{}}^{\mathrm{}}𝑑\tau ^{}`$, where the lagrangian density $``$ is given by

$`=`$ $`|\dot{P}_1|^2+|\dot{P}_2|^2i(\delta v)\left(P_1^{}\dot{P}_1P_1\dot{P}_1^{}\right)`$ (7)
$`i(\delta +v)\left(P_2^{}\dot{P}_2+P_2\dot{P}_2^{}\right)+2(q\beta )|P_1|^2`$
$`+2(q+\beta )|P_2|^2|P_1|^4|P_2|^42A|P_1|^2|P_2|^2`$
$`B\left(P_1^2P_2^{\mathrm{\hspace{0.17em}2}}+P_1^2P_2^2\right).`$
We now assume the following ansatz for $`P_{1,2}`$

$`P_i=a_isech(b_i\tau ^{})\mathrm{exp}(ic_i\tau ^{}),`$ $`i=1,2.`$ (8)
Where the variational parameters $`a_{1,2}`$, $`b_{1,2}`$, and $`c_{1,2}`$ are assumed to have real values. The above ansatz is inspired from solutions 2, and 3, as well as from the work of Paré . We know that the exact solution has a complex phase chirp  induced by the four wave mixing term (last term in Eqs.2), which we neglect for simplicity but also because from previous numerical calculations we observed that the phase dependence on $`\tau ^{}`$ was mainly lineal .

By introducing Eqs.8 into Eqs.7, and applying the Euler-Lagrange equations to the averaged Lagrangian one obtains, after some straightforward algebra, a set of six nonlinear coupled equations,viz,

$$c_1=\delta +v+Bb_1a_2^2I(c_{})$$
(9)
$$c_2=\delta +vBb_2a_1^2I(c_{})$$
(10)
$$\frac{b_1}{3}+\frac{2D_1}{b_1}\frac{4a_1^2}{3b_1}A\frac{a_2^2}{b_2}m_0(\eta )Ba_2^2M(b_1,b_2,c_{})=0$$
(11)
$$\frac{b_1}{3}\frac{2D_1}{b_1}+\frac{2a_1^2}{3b_1}A\frac{a_2^2\eta }{b_2}\dot{m}_0(\eta )Ba_2^2b_1\frac{}{b_1}M(b_1,b_2,c_{})=0$$
(12)
$$\frac{b_2}{3}+\frac{2D_2}{b_2}\frac{4a_2^2}{3b_2}A\frac{a_1^2}{b_1}m_0(\eta ^1)Ba_1^2M(b_1,b_2,c_{})=0$$
(13)
$$\frac{b_2}{3}\frac{2D_2}{b_2}+\frac{2a_2^2}{3b_2}A\frac{a_1^2}{b_1\eta }\dot{m}_0(\eta ^1)Ba_1^2b_2\frac{}{b_2}M(b_1,b_2,c_{})=0$$
(14)
where

$$I(x)=_{\mathrm{}}^{\mathrm{}}\tau \mathrm{sin}(2x\tau )sech^2(b_1\tau )sech^2(b_2\tau )𝑑\tau $$
(15)
$$c_{}=c_2c_1=2\delta \frac{b_1b_2}{2}QBI(c_{})$$
(16)
$$m(x)=_{\mathrm{}}^{\mathrm{}}sech^2tsech^2(xt)𝑑t$$
(17)
$$M(b_1,b_2,c_{})=_{\mathrm{}}^{\mathrm{}}cos(2c_{}t)sech^2(b_1t)sech^2(b_2t)𝑑t.$$
(18)
$`D_1=q\beta +c_1^2/2+(\delta v)c_1`$, $`D_2=q+\beta +c_2^2/2(\delta +v)c_2`$, $`\eta =b_1/b_2`$, and $`Q=2(a_1^2/b_1+a_2^2/b_2)`$

In order to simplify the set of equations 9 \- 14, we need to do some assumptions and approximations. We first take into account  that the main contribution to $`c_{}`$ becomes from $`2\delta `$, while the integral is a minor correction. Therefore, using the first order Börn approximation  $`c_{}`$ becomes

$`c_{}2\delta +{\displaystyle \frac{b_1b_2}{2}}QBI(2\delta )`$ (19)
Introducing Eq.19 into Eqs.11 \- 14 we get a reduced set of four nonlinear equations. The next step must consist in finding accurate analytical expressions for the integrals $`I`$, $`m`$, and $`M`$ which can not be exactly integrated. Following ref. we expand $`m`$ in a Taylor series around $`\eta =1`$ up to second order, and obtain:

$`m(x)={\displaystyle \frac{4}{3}}+{\displaystyle \frac{2}{3}}(1x)+{\displaystyle \frac{2}{3}}(1{\displaystyle \frac{\pi ^2}{15}})(1x)^2`$
$`\dot{m}(x)={\displaystyle \frac{2}{3}}{\displaystyle \frac{4}{3}}(1{\displaystyle \frac{\pi ^2}{15}})(1x)+({\displaystyle \frac{\pi ^210}{5}})(1x)^2`$ (20)
The same Taylor expansion for $`M`$ or $`I`$ gives much more complicated expressions in terms of Polygamma functions , which should be evaluated numerically. Instead of this approach we choose to approximate $`M`$ and $`I`$ by replacing $`b_i`$ by $`(b_1+b_2)/2`$. In this way the integrals can be analytically performed to give:

$`M{\displaystyle \frac{16\pi c_{}}{3b^2}}\left(1+\left({\displaystyle \frac{2c_{}}{b}}\right)^2\right)csch\left[{\displaystyle \frac{2\pi c_{}}{b}}\right]`$ (21)
$`{\displaystyle \frac{M}{b_i}}={\displaystyle \frac{108\pi c_{}^3}{3b^5}}csch\left[{\displaystyle \frac{2\pi c_{}}{b}}\right]{\displaystyle \frac{32\pi c_{}[b^2+4c_{}^2]}{3b^5}}csch\left[{\displaystyle \frac{2\pi c_{}}{b}}\right]`$
$`+{\displaystyle \frac{32\pi ^2c_{}^2[b^2+4c_{}^2]}{3b^6}}coth\left[{\displaystyle \frac{2\pi c_{}}{b}}\right]csch\left[{\displaystyle \frac{2\pi c_{}}{b}}\right]`$ (22)
$`I(x){\displaystyle \frac{16\pi x[b^2+x^2]}{3b^5}}\left[\pi coth({\displaystyle \frac{2\pi x}{b}}){\displaystyle \frac{b}{2x}}{\displaystyle \frac{4bx}{b^2+4x^2}}\right]csch[{\displaystyle \frac{2\pi x}{b}}],`$ (23)
where $`b=b_1+b_2`$.

From here the procedure goes as follow. Given the equation coefficients $`\delta `$, $`\beta `$, $`A(=2/3)`$, and $`B(=1/3)`$, together with the parameters $`v`$ and $`q`$, we solve the Eqs.11 -14 to obtain $`b_{1,2}`$, and $`a_{1,2}`$. Then $`c_{1,2}`$ are obtained from Eqs.9\- 10.

## 4 Numerical results

We have numerically solved the four coupled nonlinear algebraic equations by a Powell hybrid method . In all cases a unique solution was instantaneoulsy obtained. In order to check the validity of the approximations used to solve the integrals we sometimes evaluated numerically $`M`$, $`m`$, and $`I`$ at each step getting almost identical results but with a much higher CPU time consumption.

Fig.1 shows the dependence of the soliton energy $`Q`$ on the nonlinear propagation constant $`q`$. This diagram was obtained by means of the above described variational method. The solid lines represent the solutions given by Eqs. 4. Different sets of almost parallel curves correspond to different values of the parameter $`\delta `$ whilst the value of $`v`$ hardly influences the shape of these curves, fixing solely the minimum allowed value of $`q`$. These curves coincide totally with those shown in fig.1(a) of Ref., which were obtained exactly by numerically solving the propagation equation for certain values of $`q`$ whatever $`v`$ was.

In addition to the solutions approximated by Eq.8, where $`a_{1,2}`$ were assumed to be real, and therefore represent solutions where $`U`$ and $`V`$ are in phase, there exist solutions where $`U`$ and $`V`$ are $`\pi /2`$ out of phase . They can be obtained directly from our variational approach if we do not imposed $`a_{1,2}`$ to be real. Alternatively our approach can account for them just by changing the ansatz (Eq.8 for $`P_i`$). If we choose the following one:

$`\stackrel{~}{P}_1=a_1sech(b_1\tau ^{})exp(ic_1\tau ^{})`$
$`\stackrel{~}{P}_2=a_2sech(b_2\tau ^{})exp(i(c_2\tau ^{}+\pi /2)),`$ (24)
and repeat the process, it is easily realized that the only change produced in the nonlinear algebraic equations is to replace $`B`$ by $`B`$. This indicates that the existence of these independent solutions is intimately related to the inclusion of the four wave mixing term in the propagation equations.

Figure 2 shows the fraction of energy in the slow mode ($`Q_1/Q`$) versus the total energy ($`Q`$) of the two-parameter family of solutions for $`\delta =\beta =1`$ as obtained from our variational method. In (a) and (b) $`P_i`$ is given by Eq.8 and Eq. 4 respectively. The results can be compared with the exact ones in Fig. 1 of Ref., which were found numerically. The symmetry of these curves is a result of the symmetry that possesses the Lagrangian under:

$`q\stackrel{~}{q}=q+{\displaystyle \frac{\stackrel{~}{v}^2v^2}{2}},v\stackrel{~}{v}={\displaystyle \frac{2\beta }{\delta }}v,`$
$`c_{2,1}\stackrel{~}{c}_{2,1}={\displaystyle \frac{2\beta }{\delta }}c_{1,2},a_2a_1,b_2b_1`$ (25)
Figs. 1, and 2 illustrate the accuracy of the variational approach as developed in the section 3. We must remark that solving the set of four nonlinear algebraic equations is an easy task that can be done really fast. In our computer (Alpha Dec 2100/500) we obtain thousands of these variational solutions in just a few CPU seconds.

Fig. 3 shows the variation of the width ratio ($`\eta =b_1/b_2`$) vs. the propagation constant $`q`$ for the stationary solutions when $`U`$ and $`V`$ are (a) in phase, and (b) in quadrature. We take $`\beta =\delta =1`$, and three values of $`v`$ are considered, which are written close to the corresponding curves. In all cases $`\eta `$ tends to $`1`$ as $`q`$ increases, i.e. for high values of $`Q`$ (see Fig. 1). The departure of $`\eta `$ from unity is higher when $`U`$ and $`V`$ are in quadrature than when they are in phase. In any case it is small except around the minimum allowed value of $`q`$, where the curves emerge and one of the components ($`V`$ in the cases of figure 3) is almost zero. This figure serves also to verify a posteriori that our approximations of integrals $`I`$, $`m`$ and $`M`$ had good basis.

Similarly Fig.4 shows (a) $`c_2c_1`$ and (b) $`c_2`$ vs. $`q`$ for the same values of the parameters. As expected, at the minimum value allowed for $`q`$, where the solution is a fast soliton (Eq.3), $`c_2`$ is exactly $`\delta +v`$. As it occurred for $`Q`$ vs. $`q`$, the value of $`c_{}`$ hardly depends on $`v`$, being mainly determined by the value of $`q`$. It is also remarkable that while the solutions for $`U`$ and $`V`$ in phase decrease its frequency difference (i.e. $`c_{}`$) as we increase $`q`$ (and therefore $`Q`$), the opposite happens for the solutions with $`U`$ and $`V`$ in quadrature. In these last cases $`c_{}`$ becomes more different to $`2\delta `$ than when $`U`$ and $`V`$ are in phase. The same happened with $`\eta `$ respect to $`1`$ (see Fig. 3). In fact we observe that the variational solutions were much more accurate for the solutions with $`U`$ and $`V`$ in phase than for those in quadrature.

Fig.5 shows the variation of the widths and peak amplitudes vs. $`q`$ for the solutions with $`U`$ and $`V`$ in phase. It can be seen that for a given value of $`\delta `$ and $`\beta `$, the widths depend almost exclusively on $`q`$, being almost completely independent on $`v`$. On the other hand the asymptotic behavior of $`a_1/a_2`$ does not depend on $`q`$ but only on $`v`$ (see Fig.5c).

We have numerically solved Eqs.2, taking as initial conditions the variational solutions corresponding to different values of the parameters $`v`$ and $`q`$ and fixed values of $`\delta `$ and $`\beta `$, viz. $`\beta =1=\delta `$. In all cases that we propagate in-phase solutions, we obtained stable propagation, whilst $`\pi /2`$-dephased solutions were sometimes unstable. Fig. 6 shows the stable propagation of the in-phase variational solution for $`q=3.2`$ and $`v=0.9`$, whereas Fig.7 shows the same for the $`\pi /2`$-dephased solution corresponding to $`q=6`$ and $`v=0.9`$. Contrary to the previous case, a small difference can be appreciated between the variational solution and the stationary one. The exact stable stationary solution (which of course is $`\pi /2`$-dephased) is reached after a very short propagation distance. Finally Fig. 8 illustrates the unstable behavior of a $`\pi /2`$-dephased solution ($`q=7.2`$, $`v=0.8`$). The solution oscillates around the stationary solution with increasing amplitude emitting a small quantity of radiation. Eventually, after the emisssion of a large quantity of radiation, a stable solution is reached (not shown in the figure). In Figs. 6-8 $`\tau ^{\prime \prime }=\tau +\delta `$.

## 5 Summary

We have developed an accurate variational approach to derive analytical approximations of the coupled pulselike stationary solutions in birefringent fibers in the most general case when differences in both phase velocities and group velocities between the two components and fastly oscillatory terms are taken into account (they are important in the subpicosecond regime). As particular cases it includes those where any term in Eqs.2 can be neglected. We have shown that the difference between the central frequencies of the components, their widht and their energies are almost independent on $`v`$, being mainly determined by $`q`$. We have also shown that the ratio $`a_1/a_2`$ is mainly determined by $`v`$. The stability of these solutions has been briefly considered, by numerically propagating them. In all the cases we propagated solution in-phase, they happened to be stable whilst those $`\pi /2`$ out of phase present some intervals of stability. A global stability analysis of these solutions remains to be done. This variational method could be useful to design soliton-dragging logic gates .

Acknowledgments

This work was supported by the Comunidad de Madrid under contract 06T/039/96 and by the CICyT under contract TIC95-0563-03. N.J.R-F acknowledges a grant from the CSIC.


text:
# Electron Double Refraction in Hybrid Systems with Rashba Spin-Orbit Coupling

## I Introduction

The most popular spin-filter device, based on Rashba Effect,  has been proposed by Datta and Das  as the electronic analog of an electro-optic modulator. The idea is to design a structure in which the spin behaves as the polarization of the light. When the light traverses an electro-optic material, two perpendicular polarizations accumulate different phases shifts and, when the beam emerges into the analyzer, the two components interfere with each other. Polarizing the light at 45 in the plane $`(y,z),`$ orthogonal to the direction propagation $`x`$, the output power collected by the analyzer, that is oriented in the same way of the input filter, is given by

$$P_0=\frac{1}{4}\left|\overline{\underset{¯}{\begin{array}{cc}1\hfill & 1\hfill \end{array}}}\left|\begin{array}{c}e^{ik_1L}\hfill \\ e^{ik_2L}\hfill \end{array}\right|\right|^2=\mathrm{cos}^2\frac{\left(k_1k_2\right)L}{2}.$$
(1)
A gate voltage controls the differential phase shift $`\left(k_1k_2\right)L`$ and allows for a modulation of the output.

In the spintronic analog the role of the electro-optic material is played by a strip of 2DEG (two dimensional electron gas) in which magnetized contacts inject and collect electrons with a specific spin orientation. The feasibility of spin injection at a ferromagnet–semiconductor interface has been experimentally demonstrated.

An electrical field $``$ (in $`y`$ direction) acts on the strip and the electrons move along $`x`$ direction. We suppose that the motion happens in a nanostructure at low temperature so that the electronic phase coherence is maintained. The velocities of the charge carriers are of the order $`10^8`$m/sec or larger and a magnetic field (directed in $`z`$ direction) appears in the rest reference frame of the charges. This kind of spin-orbit effect is known as Rashba effect. The Hamiltonian spin-orbit term has the form

$$H_{SO}=\frac{g\left|e\right|}{m^2c^2}\left(\stackrel{}{p}\times \stackrel{}{}\right)\stackrel{}{s}=\frac{g\left|e\right|}{m^2c^2}\left(\stackrel{}{s}\times \stackrel{}{p}\right)\stackrel{}{},$$
(2)
where $`g`$ is the giromagnetic ratio of the spin $`\stackrel{}{s}`$ and $`m`$ stands for the effective mass of the electron charge $`\left|e\right|`$. Introducing the Pauli matrices $`\sigma _x,\sigma _y,\sigma _z`$, we get

$`H_{so}={\displaystyle \frac{\eta }{\mathrm{}}}\left(\sigma _zp_x\sigma _xp_z\right)`$

where $`\eta =\left(g\left|e\right|\mathrm{}^2/2m^2c^2\right)`$.

The motion in $`x`$ direction is considered taking momentum eigenvalues $`p_z=0`$ and $`p_x=\mathrm{}k^{}`$. In this subspace the spin component in $`z`$ direction is a motion constant and the energy eigenvalues of spin up (+) and spin down (–) states are respectively:

$$E_\pm =\frac{\mathrm{}^2}{2m}k^{\mathrm{\hspace{0.17em}\hspace{0.17em}2}}\pm \eta k^{}.$$
(3)
For the InGaAs/InAlAs heterostructure the spin-orbit parameter $`\eta `$ was estimated to be $`3.9\times 10^{12}`$ eV m. For positive values of the energy we introduce the vector $`k`$ such that $`E_\pm =\mathrm{}^2k^2/2m`$. The two spin orientations have two different values of the momentum

$$k^{}=k_\pm =\sqrt{k^2+\frac{m^2\eta ^2}{\mathrm{}^4}}\frac{m\eta }{\mathrm{}^2}.$$
(4)
At high values of the energy $`(Em\eta ^2/\mathrm{}^2)`$ the two values of $`k^{}`$ have a difference of $`2m\eta /\mathrm{}^2`$. Preparing the electron in the state $`\frac{1}{\sqrt{2}}\left|\begin{array}{c}1\hfill \\ 1\hfill \end{array}\right|`$, when it traverses a distance $`L`$, the probability that it can be found again in this state is given by:

$$P=\frac{1}{4}\left|\overline{\underset{¯}{\begin{array}{cc}1\hfill & 1\hfill \end{array}}}\left|\begin{array}{c}e^{ik_+L}\hfill \\ e^{ik_{}L}\hfill \end{array}\right|\right|^2=\mathrm{cos}^2\frac{2m\eta L}{\mathrm{}^2}.$$
(5)
The Datta and Das ideas have inspired some investigation on spintronic devices that exhibit spin-valves effects. Recently it has been shown that an oscillatory spin-filtering due to a spin-dependent conductance can not be obtained by a single interface with the Rashba coupling on one side, unlikely what previously stated. This issue stems from the boundary conditions that guarantee the continuity of the probability current density perpendicular to the interface. These conditions imply the same transmission amplitudes for spin up and spin down and so the effect is absent. The aim of this paper is to show that the double refraction arising when the electron incidence is out of normal mixes the in and out spin states allowing for an oscillatory behavior of the conductance with respect to ingoing spin polarization.

## II Spin scattering against interface

In this section we will study an hybrid system with an $`H_{so}`$ coupling described by the Hamiltonian:

$$H=\stackrel{}{p}\frac{1}{2m\left(x\right)}\stackrel{}{p}+\frac{\eta \left(x\right)}{\mathrm{}}\left(\sigma _zp_x\sigma _xp_z\right)i\sigma _z\frac{1}{2}\frac{\eta \left(x\right)}{x}+U\delta \left(x\right).$$
(6)
The spatial variation of the spin-orbit coupling $`\eta `$ and of the effective mass $`m`$ on passing from one side of the interfaces to the other are taken into account in such a way to ensure Hamiltonian hermiticity. The kinetic energy and the spin-orbit $`H_{so}`$ contain the momentum operator and have been symmetrized in Eq.(6).

We assume that the parameters are piecewise constant

$`m\left(x\right)`$ $`=`$ $`m_F\vartheta \left(x\right)+m_S\vartheta \left(x\right)`$ (7)
$`\eta \left(x\right)`$ $`=`$ $`\eta \vartheta \left(x\right),`$ (8)
where $`\vartheta \left(x\right)`$ is the step function. We have added a term $`U\delta \left(x\right)`$ to control the transparency of the interface. The spinor eigenstate of $`H,\psi `$ is continuous while its derivative has a discontinuity fixed by the strength of the Dirac delta in $`x=0`$:

$`\psi \left(0+\right)`$ $`=`$ $`\psi \left(0\right)`$ (9)
$`{\displaystyle \frac{\psi \left(0+\right)}{x}}\mu {\displaystyle \frac{\psi \left(0\right)}{x}}`$ $`=`$ $`\left(uik_0\sigma _z\right)\psi \left(0\right)`$ (10)
where $`\mu =m_S/m_F`$, $`u=2mU/\mathrm{}^2`$ and $`k_0=m\eta /\mathrm{}^2.`$

The free 2DEG with Rashba term occupying the whole $`xz`$ plane has the spinors:

$`\psi _+`$ $`=`$ $`\mathrm{exp}i\left(k_1x+k_2z\right)\left(\begin{array}{c}\mathrm{sin}\theta \hfill \\ \mathrm{cos}\theta \hfill \end{array}\right)`$ (13)
$`\psi _{}`$ $`=`$ $`\mathrm{exp}i\left(k_1x+k_2z\right)\left(\begin{array}{c}\mathrm{cos}\theta \hfill \\ \mathrm{sin}\theta \hfill \end{array}\right)`$ (16)
as energy eigenstates of the eigenvalues:

$$E_\pm =\frac{\mathrm{}^2}{2m_S}\left(k_1^2+k_2^2\right)\pm \eta \sqrt{k_1^2+k_2^2},$$
(17)
where $`\mathrm{}k_1`$ and $`\mathrm{}k_2`$ are, respectively, the eigenvalues of the momentum components $`p_x`$ and $`p_z`$, and

$$\theta =\mathrm{arctan}\left[\frac{k_1}{k_2}+\sqrt{\frac{k_1^2}{k_2^2}+1}\right].$$
(18)
If $`k^{}=\sqrt{k_1^2+k_2^2}`$ is the modulus of the momentum and $`\varphi `$ gives the direction of the in-plane motion ( $`k_1/k_2=\mathrm{cot}\varphi `$ ), then

$$\theta =\frac{\pi }{2}\frac{\varphi }{2}$$
(19)
and

$$E_\pm =\frac{\mathrm{}^2}{2m_S}\left(k^{\mathrm{\hspace{0.17em}\hspace{0.17em}2}}\pm 2k_0k^{}\right).$$
(20)
We note that the two spinors $`\psi _+`$ and $`\psi _{}`$ are one the time reversed of the other. The time reversal operator $`\widehat{T}`$

$`\widehat{T}\left(\begin{array}{c}\psi _1\hfill \\ \psi _2\hfill \end{array}\right)=\left(\begin{array}{c}\psi _2^{}\hfill \\ \psi _1^{}\hfill \end{array}\right)`$

commutes with $`H_{SO}`$: $`[\widehat{T},H_{SO}]=0`$. The degeneracy with respect to $`k^{}`$ is lifted but the Rashba term is not able to produce a spontaneous spin polarization of the electron states: any given energy value fixes two different values of the modulus $`k^{}`$, leaving undetermined the spin polarization. However we emphasize that the direction of the wave vector $`\stackrel{}{k}^{}`$ fixes the spin polarization as the equation (19) shows. If we choose the in-plane motion direction then we put the electron in a definite spin polarization state. If $`\stackrel{}{k}^{}`$ is along $`x`$ direction then $`\varphi =0`$ and $`\psi _+,\psi _{}`$ are the up and down spins in $`z`$ direction. Including the spin-orbit interaction in the Hamiltonian the double group is the new space group  and a space rotation of $`4\pi `$ is needed to achieve the same spinor.

Let be

$`E={\displaystyle \frac{\mathrm{}^2}{2m_F}}k^2`$

the electron energy in the ferromagnet. When the electron goes into the 2DEG region its wave vector $`k^{}`$ becomes

$$k^{}=\sqrt{\mu k^2+k_0^2}k_0=k_\pm ,$$
(21)
where the index $`\pm `$ refers to the two branches with the same energy $`E`$

$`E_\pm ={\displaystyle \frac{\mathrm{}^2}{2m_S}}\left(k_\pm ^2\pm 2k_0k_\pm \right).`$

The hybrid system is invariant for translation along $`z`$ direction and the component $`p_z`$ of the momentum is conserved. If $`\alpha `$ is the angle of $`\stackrel{}{k}_+`$ with the $`x`$ axis, $`\beta `$ and $`\gamma `$ the corresponding angles of $`\stackrel{}{k}_{}`$ and $`\stackrel{}{k}`$, respectively, the momentum conservation implies that

$$k_+\mathrm{sin}\alpha =k_{}\mathrm{sin}\beta =k\mathrm{sin}\gamma .$$
(22)
The figure 1 shows the output angles $`\alpha `$ and $`\beta `$. Only when the incidence is normal, i.e. $`\gamma =0`$, the electron moves in the 2DEG in the same direction with $`\alpha =\beta =0`$. In all other cases, i.e. $`\gamma >0`$, the outgoing states $`+`$ and $``$ go along two different directions. This phenomenon is analogous to the double refraction that appears in biaxial crystals. Again the spin of the electrons behaves as the polarization of the light. We remember that the birefringence arises when the characteristics of electromagnetic propagation depend on the directions of propagation and polarization of the wave.

The mode $`+`$ has the limit angle

$$\gamma _0=\mathrm{arcsin}\frac{k_+}{k},$$
(23)
so that for $`\gamma >\gamma _0`$ this mode is totally reflected: it vanishes exponentially for $`x>0`$. Here and in the following we take $`0<\mu <1`$ because the effective mass in the 2DEG is less than in the ferromagnet. When $`k/k_0<2/\left(1\mu \right)`$ the mode $``$ is always transmitted up to grazing incidence at $`\gamma =\pi /2`$. Increasing the kinetic energy with respect to spin-orbit coupling, when $`k/k_0>2/\left(1\mu \right)`$, a second limit angle appears:

$$\gamma _1=\mathrm{arcsin}\frac{k_{}}{k}>\gamma _0$$
(24)
and, for $`\gamma >\gamma _1`$, we have the total reflection (both the modes vanish for $`x>0`$). When the strength of spin-orbit coupling goes to zero, $`\gamma _0`$ and $`\gamma _1`$ tend to the common limit $`\mathrm{arcsin}\sqrt{\mu }`$.

Lighter the effective mass within the 2DEG is, nearer to normal incidence the propagation directions $`\alpha `$ and $`\beta `$ into Rashba region $`x>0`$ are. The figure 2 shows the limit angles as a function of $`k/k_0`$.

The incoming spinor

$$\psi _i=\mathrm{exp}\left(ik\left(x\mathrm{cos}\gamma +z\mathrm{sin}\gamma \right)\right)\left(\begin{array}{c}\mathrm{cos}\delta \hfill \\ \mathrm{sin}\delta \hfill \end{array}\right)$$
(25)
is reflected at the interface $`x=0`$ as

$$\psi _r=\mathrm{exp}\left(ik\left(x\mathrm{cos}\gamma +z\mathrm{sin}\gamma \right)\right)\left(\begin{array}{c}r_+\hfill \\ r_{}\hfill \end{array}\right)$$
(26)
and transmitted at $`x>0`$ in both the modes $`+`$ and $``$ as

$`\psi _t`$ $`=`$ $`t_+\mathrm{exp}\left(ik_+\left(x\mathrm{cos}\alpha +z\mathrm{sin}\alpha \right)\right)\left(\begin{array}{c}\mathrm{cos}\alpha /2\\ \mathrm{sin}\alpha /2\end{array}\right)+`$ (32)
$`\dot{t}_{}\mathrm{exp}\left(ik_{}\left(x\mathrm{cos}\beta +z\mathrm{sin}\beta \right)\right)\left(\begin{array}{c}\mathrm{sin}\beta /2\\ \mathrm{cos}\beta /2\end{array}\right),`$
where $`\delta `$ fixes the spin polarization of electron within the ferromagnet. When $`\gamma =0`$ a spin up (along $`z`$ direction) goes in the mode $`+,`$ while the spin down propagates in the mode $``$ at $`x>0.`$ In this case Zülike et al. and Molenkamp et al.  have shown that $`t_+=t_{}`$ and the interface is not able to filter the spin. We note that out of the normal incidence with $`\gamma >0`$ the scattering changes the spin polarization. The transmitted amplitudes $`t_+,t_{}`$ and the reflected ones $`r_+,r_{}`$ are determined by the boundary conditions (9) as functions of $`k,k_0,u,\delta `$ and $`\gamma `$. Solving the system

$`t_+\mathrm{cos}{\displaystyle \frac{\alpha }{2}}t_{}\mathrm{sin}{\displaystyle \frac{\beta }{2}}=\mathrm{cos}\delta +r_+`$

$`t_+\mathrm{sin}{\displaystyle \frac{\alpha }{2}}+t_{}\mathrm{cos}{\displaystyle \frac{\beta }{2}}=\mathrm{sin}\delta +r_{}`$

$`k_+t_+\mathrm{cos}\alpha \mathrm{cos}{\displaystyle \frac{\alpha }{2}}k_{}t_{}\mathrm{cos}\beta \mathrm{sin}{\displaystyle \frac{\beta }{2}}\mu k\mathrm{cos}\gamma \left(\mathrm{cos}\delta r_+\right)=`$

$$\left(k_0+iu\right)\left(\mathrm{cos}\delta +r_+\right)$$
(33)
$`k_+t_+\mathrm{cos}\alpha \mathrm{sin}{\displaystyle \frac{\alpha }{2}}+k_{}t_{}\mathrm{cos}\beta \mathrm{cos}{\displaystyle \frac{\beta }{2}}\mu k\mathrm{cos}\gamma \left(\mathrm{sin}\delta r_{}\right)=`$

$`\left(k_0iu\right)\left(\mathrm{sin}\delta +r_{}\right).`$

We get

$`r_+=(C_+A_{}\mathrm{cos}\alpha /2C_{}A_+\mathrm{sin}\beta /2)/D`$

$`r_{}=(C_{}A_{++}\mathrm{cos}\beta /2C_+A_+\mathrm{sin}\alpha /2)/D`$

$$t_+=\left[\left(\mathrm{cos}\delta +r_+\right)\mathrm{cos}\beta /2+\left(\mathrm{sin}\delta +r_{}\right)\mathrm{sin}\beta /2\right]/\mathrm{cos}\frac{\alpha \beta }{2}$$
(34)
$`t_{}=\left[\left(\mathrm{sin}\delta +r_{}\right)\mathrm{cos}\alpha /2\left(\mathrm{cos}\delta +r_+\right)\mathrm{sin}\alpha /2\right]/\mathrm{cos}{\displaystyle \frac{\alpha \beta }{2}}`$

with

$`A_{++}`$ $`=`$ $`k_+\mathrm{cos}\alpha +\mu k\mathrm{cos}\gamma +k_0+iu`$
$`A_+`$ $`=`$ $`k_+\mathrm{cos}\alpha +\mu k\mathrm{cos}\gamma k_0+iu`$
$`A_+`$ $`=`$ $`k_{}\mathrm{cos}\beta \mu k\mathrm{cos}\gamma k_0iu`$
$`A_{}`$ $`=`$ $`k_{}\mathrm{cos}\beta +\mu k\mathrm{cos}\gamma k_0+iu`$
$`C_+`$ $`=`$ $`\left(k_+\mathrm{cos}\alpha +\mu k\mathrm{cos}\gamma k_0iu\right)\mathrm{cos}\delta \mathrm{cos}{\displaystyle \frac{\beta }{2}}+`$
$`\left(k_+\mathrm{cos}\alpha +\mu k\mathrm{cos}\gamma +k_0iu\right)\mathrm{sin}\delta \mathrm{sin}{\displaystyle \frac{\beta }{2}}`$
$`C_{}`$ $`=`$ $`\left(k_{}\mathrm{cos}\beta \mu k\mathrm{cos}\gamma +k_0+iu\right)\mathrm{cos}\delta \mathrm{sin}{\displaystyle \frac{\alpha }{2}}+`$
$`\left(k_{}\mathrm{cos}\beta +\mu k\mathrm{cos}\gamma +k_0iu\right)\mathrm{sin}\delta \mathrm{cos}{\displaystyle \frac{\alpha }{2}}`$
$`D=A_{++}A_{}\mathrm{cos}{\displaystyle \frac{\beta }{2}}\mathrm{cos}{\displaystyle \frac{\alpha }{2}}A_+A_+\mathrm{sin}{\displaystyle \frac{\beta }{2}}\mathrm{sin}{\displaystyle \frac{\alpha }{2}}.`$

We note that when $`\gamma >\gamma _0`$ then

$`\mathrm{sin}\alpha ={\displaystyle \frac{k}{k_+}}\mathrm{sin}\gamma >1,`$

whose solution in $`\alpha `$ is

$`\alpha ={\displaystyle \frac{\pi }{2}}+i\alpha ^{};\mathrm{sin}\alpha =\mathrm{cosh}\alpha ^{};\mathrm{cos}\alpha =i\mathrm{sinh}\alpha ^{}.`$

The mode $`+`$ becomes a vanishing wave along $`x`$ axis:

$`t_+\mathrm{exp}\left(k_+x\mathrm{sinh}\alpha ^{}\right)\mathrm{exp}\left(ik_+z\mathrm{cosh}\alpha ^{}\right)\left(\begin{array}{c}\mathrm{cos}\left(\pi /4+i\alpha ^{}/2\right)\hfill \\ \mathrm{sin}\left(\pi /4+i\alpha ^{}/2\right)\hfill \end{array}\right).`$

When $`\gamma >\gamma _1,`$ $`\beta =\pi /2+i\beta ^{}`$ and both the modes are damped within the 2DEG: the incident wave is totally reflected.

At normal incidence

$`\gamma =\alpha =\beta =0`$

and

$`t_+`$ $`=`$ $`{\displaystyle \frac{2\mu k\mathrm{cos}\delta }{k_++k_0+iu+\mu k}}`$
$`t_{}`$ $`=`$ $`{\displaystyle \frac{2\mu k\mathrm{sin}\delta }{k_.k_0+iu+\mu k}}.`$
Since

$$k_++k_0=k_.k_0=\sqrt{\mu k^2+k_0^2}$$
(35)
the transmitted spinor is:

$`\psi _t`$ $`=`$ $`{\displaystyle \frac{2\mu k}{\sqrt{\mu k^2+k_0^2}iu+\mu k}}(e^{ik_+x}\left(\begin{array}{c}\mathrm{cos}\delta \hfill \\ 0\hfill \end{array}\right)`$
$`+e^{ik_{}x}\left(\begin{array}{c}0\hfill \\ \mathrm{sin}\delta \hfill \end{array}\right))`$
and the interference between the modes $`+`$ and $``$ at the interface in $`x=0`$ is lost. If $`\delta =\pi /4`$, projecting $`\psi _t`$ on the input spinor, we get

$`\left|\psi _t^{}\times \left(\begin{array}{c}1/\sqrt{2}\hfill \\ 1/\sqrt{2}\hfill \end{array}\right)\right|^2\mathrm{cos}^2{\displaystyle \frac{k_+k_{}}{2}}x`$

that is the Datta and Das  modulation factor. However an analyzer, i.e. a second interface, is needed to have a spin dependent transmission.

The square moduli of the transmitted amplitudes $`\left|t_\pm \left(\delta \right)\right|^2`$ are shown in fig.3 when $`\gamma `$ is between $`0`$ and $`\pi /2`$. We see how $`\left|t_\pm \left(0\right)\right|^2`$ and $`\left|t_\pm \left(\pi /2\right)\right|^2`$, and $`\left|t_\pm \left(\pi /4\right)\right|^2`$ and $`\left|t_\pm \left(3\pi /4\right)\right|^2`$ too, start from the same value for $`\gamma =0`$ but are different when the incidence angle increases towards $`\pi /2`$. The derivatives of $`\left|t_\pm \left(\delta \right)\right|^2`$ jumps at $`\gamma _0`$ and then at $`\gamma _1`$, when the character of the mode propagation changes. The traversing of the interface changes the spin polarization when $`\gamma >0`$.

We get the transmission coefficient $`T`$ from the probability current density

$$\begin{array}{cc}\stackrel{}{j}=\mathrm{}\left\{\psi ^{}\stackrel{}{p}\psi \right\}\hfill & ;x<0\hfill \\ \stackrel{}{j}=\mathrm{}\left\{\psi ^{}\left(\stackrel{}{p}+\mathrm{}k_0\widehat{y}\times \stackrel{}{\sigma }\right)\psi \right\}\hfill & ;x>0\hfill \end{array}$$
(38)
whose $`x`$component is

$`j_{xr}`$ $`=`$ $`\mathrm{}k\mathrm{cos}\gamma \left(1\left|r_+\right|^2\left|r_{}\right|^2\right)/m_F`$ (39)
$`\text{for}x`$ $`<`$ $`0\text{and}`$ (40)
$`j_{xl}`$ $`=`$ $`[\mathrm{}(k_++k_0)\mathrm{cos}\alpha |t_+|^2+`$ (42)
$`(k_{}k_0)\mathrm{cos}\beta \left|t_{}|^2\right]/m_S`$
$`\text{for }x`$ $`>`$ $`0.`$ (43)
The boundary conditions (9) assure the continuity of $`j_x`$ as can be verified by a straightforward calculation of eqs.(39). When $`\gamma <\gamma _0`$ both the modes propagate in $`x>0,`$ the only $``$ mode remains when $`\gamma _0<\gamma <\gamma _1`$.

The transmission coefficient is the ratio of $`j_{xr}`$ with the incident flux $`j_i=\mathrm{}k\mathrm{cos}\gamma /m_F`$, $`T=j_{xr}/j_i`$, while the reflection coefficient is $`R=\left(j_ij_{xr}\right)/j_i`$:

$`T(\delta ,\gamma )`$ $`=`$ $`((k_++k_0)\mathrm{cos}\alpha |t_+|^2\vartheta (\gamma _0\gamma )+`$ (45)
$`(k_{}k_0)\mathrm{cos}\beta \left|t_{}|^2\vartheta (\gamma _1\gamma )\right)/\mu k\mathrm{cos}\gamma `$
$`R(\delta ,\gamma )`$ $`=`$ $`\left|r_+\right|^2+\left|r_{}\right|^2`$ (46)
and when $`\gamma `$ overcomes $`\gamma _1`$ , $`T(\delta ,\gamma )=0`$ and $`R(\delta ,\gamma )=1`$. The flux is conserved because in all the cases

$`T(\delta ,\gamma )+R(\delta ,\gamma )=1.`$

The transmission coefficient as a function of $`\gamma `$ has a first higher step up to $`\gamma _0`$ followed by a lower step that ends in $`\gamma _1`$. The fig.4 shows how the shapes and the heights of the two steps vary with the spin polarization angle $`\delta `$. At low values of $`\mu `$, that is the electrons in 2DEG are light, the propagation in the $`x>0`$ region happens at angles nearer to the normal incidence. At equal masses $`(\mu =1)`$ the passage is allowed up to grazing incidence and the steps appear more squared. We note that the second step tends to disappear around $`\delta =\pi /4`$ and has the maximum height around $`\delta =3\pi /4.`$ The fig.5 refers to the case of an higher Fermi wave vector $`k.`$ Obviously when $`k/k_0\mathrm{}`$, $`T=1`$ for $`\gamma `$ from $`0`$ to $`\pi /2`$ but the second step is again visible for $`k`$ greater then $`k_0`$ of two magnitude orders.

## III Quantum Point Contact Conductance

The previously described double refraction can affect the conductance of a ballistic quantum point contact.

Let a constriction of width $`W`$ separate the ferromagnet and the 2DEG that behave as two perfect reservoirs at the Fermi energy:

$`E_F={\displaystyle \frac{\mathrm{}^2k^2}{m_F}}=E_\pm ={\displaystyle \frac{\mathrm{}^2}{2m_S}}\left(k^{\mathrm{\hspace{0.17em}\hspace{0.17em}2}}\pm 2k_0k^{}\right).`$

The electron motion within the hybrid system is assumed to be ballistic; that is the electronic mean free path is much longer than the size $`W`$ of the point contact. The Landauer-Büttiker formalism applies .The conductance $`G`$ at zero temperature is given by

$$G=\frac{e^2}{h}\underset{i}{}T_i,$$
(47)
where $`T_i`$ are the transmission coefficients for all the open channels $`i`$between the two reservoirs at the energy $`E_F`$ .In our case the index $`i`$ represents the incidence angle $`\gamma .`$

A sketch of the point contact can be found in fig.6a). The 2D Fermi circle in k-space appears in fig.6b) and only the states on its edge can carry current at zero temperature. As we have shown before, the current is transported through the point contact by the states belonging to the arch from $`\gamma _1`$ to $`\gamma _1`$ on the Fermi circle.

Quantum mechanically, the current through the point contact is equipartitioned among the $`1D`$ sub-bands, or transverse modes, in the constriction. The gap along $`k_z`$ axis between two consecutive sub-bands can be estimated of the order of $`\pi /W`$ (this is exactly the result for a square well lateral confining potential of width $`W`$). The number of states contained in the element of arch $`d\gamma `$ is then $`kd\gamma /\left(\pi /W\right).`$ The equation (47) implies that hybrid system conductance $`G`$ is

$$G=\frac{e^2}{h}_{\gamma _1}^{\gamma _1}T(\delta ,\gamma )\frac{kWd\gamma }{\pi }=\frac{e^2kW}{h}𝒢\left(\delta \right)$$
(48)
with

$$𝒢\left(\delta \right)=\frac{1}{\pi }_{\gamma _1}^{\gamma _1}T(\delta ,\gamma )𝑑\gamma .$$
(49)
An exhaustive discussion about this approach can be found in references 15 and 16. We note that the restriction to the normal incidence $`\gamma =0`$ gives

$`𝒢\left(\delta \right)={\displaystyle \frac{T(\delta ,0)}{\pi }}`$

that is the Sharvin resistance formula  used by Grundler  but that is independent on the spin polarization angle $`\delta `$.

The fig.7 shows $`𝒢\left(\delta \right)`$ for $`\delta `$ between 0 and $`\pi `$. The oscillatory behavior of the conductance allows the spin filtering with a single interface. This effect is a direct consequence of the double refraction at the interface that changes the spin state when the electron enters the region where the spin-orbit Rashba coupling works. At normal incidence the electron pass into 2DEG conserving the spin state. When a lateral confining potential is imposed to the electron gas the Q1DEG has sub-bands for which the free electron property (35) is no more valid, although the time reversal symmetry leaves the degeneracy of states with opposite value of $`k_1`$. The case of a parabolic confining potential has been studied by Governale et al.  that estimate the deformation of sub-bands and the lateral spin density. The ballistic spin-transport properties of a quasi-one-dimensional wire with a spin-orbit Rashba interaction in a finite piece of it have been studied with a numerical tight binding model by Mirales et al. They find a spin-conductance modulation.

An estimation of the strength of Rashba interaction on the conductance is given by the ratio:

$`\mathrm{\Delta }𝒢/𝒢={\displaystyle \frac{𝒢\left(3\pi /4\right)𝒢\left(\pi /4\right)}{𝒢\left(0\right)}}`$

reported in Fig.8. That is roughly the maximum relative variation of the conductance against $`k/k_0`$. We note that for $`k/k_0=100`$, $`\mathrm{\Delta }𝒢/𝒢`$ is of the order of ten per cent. We think that such a value could be detected experimentally in a quantum point contact.

## IV Conclusions

In this paper transmission across a ferromagnet/2DEG has been studied. An electron in a definite state of spin undergoes a double refraction traversing the interface analogously to what happens to the polarized light impinging the surface of a biaxial crystal. We have shown that the correct boundary conditions give rise to a spin-dependent transmission coefficient and that the normal incidence is a special case for which the dependence on spin is lost. The spin filtering occurs when the electron hits the interface in a direction out of the normal. The conductance of a point contact at the interface in ballistic transport regime within the Landauer-Büttiker formalism has been calculated. We have shown that the conductance has an oscillatory behavior with the polarization angle of the spin.

We gratefully acknowledge M.Governale for helpful suggestions about the correct way to impose the boundary conditions at the interface.

Figure Captions

1. The vectors $`\stackrel{}{k}_+,`$ $`\stackrel{}{k}_{},`$ $`\stackrel{}{k}`$ in $`k`$space and the angles $`\alpha ,\beta `$ and $`\gamma `$ that they form with $`x`$ direction normal to the interface. The two circles are the lines at the constant energy $`\mathrm{}^2k^2/2m_F.`$
2. The limit angles $`\gamma _0`$ of $`+`$ mode (dashed line) and $`\gamma _1`$ of $``$ mode (full line) for the indicated values of mass ratio $`\mu `$ as functions of $`k/k_0`$. For $`\gamma `$ above $`\gamma _1`$ the total reflection occurs.
3. The squared moduli of the transmitted amplitudes for two couples of orthogonal spin polarizations. The cusps sign the passage through the limit angles.
4. The two steps of the transmission coefficient $`T`$. The second step tends to disappear for $`\delta =\pi /4`$ and to have the same height of the first when $`\delta =3\pi /4`$.
5. The steps of $`T`$ at the higher value of $`k/k_0=100.`$ In the limit $`k/k_0\mathrm{},`$ $`T=1`$ for $`0<\gamma <\pi /2.`$
6. a)The sketch of the point contact
b) The Fermi circle in $`k`$space. The thick arch indicates the states that carry current into the point contact.
7. The conductance $`𝒢`$ as a function of the polarization spin angle $`\delta `$.
8. The relative variation of the conductance $`\mathrm{\Delta }𝒢/𝒢`$ against $`k/k_0`$.


text:
# Microscopic nonequilibrium dynamics of an inhomogeneous Bose gas beyond the Born approximation

## I Introduction

In the mean-field approximation, a Bose-condensed phase is well described by the Gross-Pitaevskii (GP) equation dalfovo . Examples of collective effects that arise at the mean-field level include the formation of vortex states murray ; haljan ; raman ; ketterle ; dalibard  and collective excitations jin96 ; edwards ; stringari . The GP equation is often sufficient to describe the dynamics at $`T=0`$. However, at finite temperatures, $`T`$ near $`T_c`$ \[the Bose-Einstein condensation (BEC) temperature\], effects due to the presence of a thermal component may be important. The hot thermal atoms interact with each other and with the condensed atoms via binary collisions. In fact, collisional dynamics is the essential mechanism for evaporative cooling leading to BEC. Finite temperature effects generate phenomena such as phase diffusion and damping of collective excitations. In order to describe these effects, inclusion of collisions arising from the thermal component becomes necessary. A generalized GP equation, for example the Hartree-Fock-Bogoliubov approach, includes collisions only as first-order energy shifts and is therefore valid only at low temperature. In the other limit, $`T>T_c`$ (high temperature), the condensate component may be negligible and the dynamics is then dominated by the thermal component. This limit is well described by the quantum-Boltzmann (QB) equation. A second order kinetic theory, which reduces to the GP and the QB equations in the limits $`T=0`$ and $`T>T_c`$, respectively, is necessary to provide a complete description.

Several different versions of quantum kinetic theories have been formulated to describe a quantum degenerate Bose system gardiner971 ; gardiner973 ; milena ; walser ; gardiner975 ; stoof ; morgan ; burnett . However certain questions remain unanswered. An important issue is the question of Markovian versus non-Markovian dynamics in a non-homogeneous system. Previous attempts to address this issue morozov ; morozov2 ; bonitz ; kremp ; bonitzbook , have been limited to the discussion of system with translational symmetry.

In the case of a sufficiently dilute nonhomogeneous system, one may consider the neglection of memory effects under the principal of the rapid attenuation of correlations. Such an ansatz will result in a Born-Markov approximation. While in this limit kinetic equations are vastly simplified, they do not conserve energy. Also, the Born-Markov theory involves a regularizing factor $`\eta 0`$ zubarev ; peletminskii . Such a limiting procedure is computationally nontrivial to implement numerically if one is dealing with a finite system with discrete energy levels. If the resolution is instead to include memory effects and initial correlations, then the full non-Markovian theory which arises is described by complex kinetic equations.

In this article we wish to address some of these issues by generalizing the kinetic theory in Ref. walser  to include effects beyond the Born approximation such as quasi-particle damping and effects arising from the finite duration of a collision. We hypothesize the form of the damping function by looking at the corresponding retarded Green’s function haug ; semkat . We show that such a perturbation theory conserves energy up to second order in the interaction strength.

We apply this theory to a simple model of a Bose gas confined in a spherical box above $`T_c`$. The simplicity of this model allows us to study the real-time (nonequilibrium) evolution of the system from some initial state to a final state of equilibrium. The equilibrium state is essentially a self-consistent steady-state solution of the kinetic equation. We perform a linear response calculation in order to verify the stability and to study the damping rates and frequencies of the steady-state solution. We show the existence of different time scales corresponding to damping and equilibration.

This article is organized as follows. In Sec. II we derive a non-Markovian generalization to the kinetic theory of Ref. walser  using the prescription of a nonequilibrium statistical operator method zubarev . Both quasi-particle damping and damping arising from the finite time of collision events will be discussed in Sec. III. This has implications on the underlying symmetries of the theory and their associated conserved quantities. In Sec. IV, we apply this generalized kinetic theory to a simple model of a dilute Bose gas confined in a spherical box and obtain a self consistent steady-state solution to the second-order kinetic theory. Finally in Sec. V, we study the response of the system to a small perturbation. This will help us to predict frequencies and damping rates of collective modes.

## II Kinetic equations

We start with the many-body Hamiltonian for a weakly interacting Bose gas given by

$$\widehat{H}=\widehat{H}^{(0)}+\widehat{H}^{(1)},$$
(1)
where $`\widehat{H}^{(0)}`$ is the single particle Hamiltonian that is defined as

$$\widehat{H}^{(0)}=\epsilon ^{12}\widehat{a}_1^{}\widehat{a}_2^{}$$
(2)
using the implicit summation convention for repeated indices. The two-body energy $`\widehat{H}^{(1)}`$ is given by

$$\widehat{H}^{(1)}=\varphi ^{1234}\widehat{a}_1^{}\widehat{a}_2^{}\widehat{a}_3^{}\widehat{a}_4^{}.$$
(3)
The bosonic operators $`\widehat{a}_1`$ and $`\widehat{a}_1^{}`$ annihilate and create a particle in a single particle state $`|1`$, respectively. The abbreviated notation $`|1`$ represents a state specified by a complete set of quantum numbers for both the motional and electronic degrees of freedom.

We assume that the particles are confined by an external trapping potential $`V_{\text{ext}}`$. Thus, the matrix elements of the single particle Hamiltonian are given by

$`\epsilon ^{12}`$ $`=`$ $`1|{\displaystyle \frac{\widehat{𝐩}^2}{2m}}+V_{\mathrm{ext}}(\widehat{𝐱})|2,`$ (4)
where $`m`$ is the mass. The binary interaction is mediated by a short range repulsive potential $`V_{\text{bin}}`$. This gives the symmetrized $`(𝒮)`$ matrix elements

$`\varphi ^{1234}`$ $`=`$ $`(𝒮)1|2|V_{\text{bin}}(\widehat{𝐱}𝟙𝟙\widehat{𝐱})|3|4.`$ (5)
In the low energy limit $`V_{\mathrm{bin}}`$ can be approximated by a contact potential with the matrix elements given by

$`\varphi ^{1234}`$ $``$ $`{\displaystyle \frac{V_0}{2}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}1|𝐱2|𝐱𝐱|3𝐱|4d^3x.`$ (6)
The interaction strength $`V_0`$, is related to the scattering length $`a_\text{S}`$ by $`V_0=4\pi \mathrm{}^2a_\text{S}/m`$.

Here, we use the well-known nonequilibrium statistical operator method zubarev ; peletminskii  to obtain an explicitly non-Markovian version of the kinetic theory walser . In this approach, the nonequilibrium state of a weakly interacting quantum gas is specified by a set of single time master variables. For our system the most important master variable is the single particle density matrix $`f(t)`$,

$$f_{12}(t)=\widehat{a}_2^{}\widehat{a}_1^{}=\mathrm{Tr}\left\{\widehat{a}_2^{}\widehat{a}_1^{}\sigma (t)\right\},$$
(7)
where $`\mathrm{}=\text{Tr}\{\mathrm{}\sigma (t)\}`$ and $`\sigma (t)`$ is the statistical many-body density operator. We focus our studies on the temperature regime above and in close proximity to the critical temperature for BEC. For this reason, we do not consider either symmetry breaking fields, $`\widehat{a}_1`$, or the anomalous fluctuations, $`\widehat{a}_1\widehat{a}_2`$. We therefore define $`\{\widehat{\gamma }_0^{}=𝟙,\widehat{\gamma }_𝕜^{}=\widehat{𝕒}_\mathrm{𝟚}^{}\widehat{𝕒}_\mathrm{𝟙}^{}|𝕜\{(\mathrm{𝟙},\mathrm{𝟚})\}\}`$ as our complete set of relevant operators. The expectation values of these operators, $`\gamma _k(t)=\widehat{\gamma }_k`$, are the only quantities that will appear in the final kinetic equations.

The time evolution of the nonequilibrium statistical operator $`\sigma (t)`$ is described by the Liouville equation with an extra source term on the right hand side:

$$\frac{d}{dt}\sigma (t)+\frac{i}{\mathrm{}}[\widehat{H},\sigma (t)]=\eta \left(\sigma (t)\sigma ^{(0)}(t)\right).$$
(8)
Such a source term breaks the time reversal symmetry of the Liouville equation and represents a convenient way to incorporate the irreversible character of macroscopic processes. The relevant distribution $`\sigma ^{(0)}(t)`$ given by

$$\sigma ^{(0)}(t)=\sigma _{\{\gamma (t)\}}^{(0)}=\mathrm{exp}\left\{\widehat{\gamma }_k\mathrm{{\rm Y}}^k(t)\right\},$$
(9)
where $`\mathrm{{\rm Y}}^k(t)`$ are the Lagrange multipliers, represents a special solution that maximizes the information entropy $`S^{}=\mathrm{Tr}\{\sigma ^{}\mathrm{log}(\sigma ^{})\}`$ for the given averages $`\gamma _k(t)`$. Furthermore, at some initial instance $`t=t_0`$ in the remote past, we can assume that $`\sigma (t_0)`$ corresponds to its noninteracting value and therefore

$$\sigma (t_0)=\sigma ^{(0)}(t_0).$$
(10)
The Lagrange multipliers $`\mathrm{{\rm Y}}^k(t)`$ are calculated from the self-consistency condition

$`\gamma _k(t)`$ $`=`$ $`\text{Tr}\{\widehat{\gamma }_k\sigma ^{(0)}(t)\}=\text{Tr}\{\widehat{\gamma }_k\sigma (t)\}.`$ (11)
This essentially enforces the Chapman-Enskog condition chapman  for the restricted set of relevant operators at all times.

From the Liouville equation (8) one can easily establish the basic equations of motion for the average values,

$`{\displaystyle \frac{d}{dt}}\gamma _k(t)`$ $`=`$ $`{\displaystyle \frac{i}{\mathrm{}}}\mathrm{Tr}\left\{[\widehat{H}^{(0)},\widehat{\gamma }_k]\sigma (t)\right\}`$ (12)
$`+`$ $`{\displaystyle \frac{i}{\mathrm{}}}\mathrm{Tr}\left\{[\widehat{H}^{(1)},\widehat{\gamma }_k]\sigma (t)\right\}.`$
The form of $`\widehat{H}^{(0)}`$ enables us to express the first trace on the right hand side of the above equation in terms of the averages $`\gamma _k(t)`$. The second trace plays the role of the “collision” term, evaluation of which requires us to seek an integral solution of the Liouville equation. But, before we proceed, it is instructive to repartition the total Hamiltonian Eq. (1) into single particle and two-particle contributions,

$`\widehat{H}`$ $`=`$ $`\widehat{\overline{H}}^{(0)}(t)+\widehat{\overline{H}}^{(1)}(t)`$ (13)
$`=`$ $`\left[\widehat{H}^{(0)}+\widehat{Q}(t)\right]+\left[\widehat{H}^{(1)}\widehat{Q}(t)\right].`$
This modification anticipates self-energy shifts

$$\widehat{Q}(t)=Q^{12}(t)\widehat{a}_1^{}\widehat{a}_2^{},$$
(14)
which will inevitably arise in the course of the calculation.

An integral solution for $`\sigma (t)`$ can then be obtained easily from the Liouville equation (8), by using the single particle time evolution operator $`\widehat{U}^{(0)}(t,t_1)`$, the boundary condition Eq. (9), and an additional partial integration. Thus, one finds

$`\sigma (t)`$ $`=`$ $`\sigma ^{(0)}(t){\displaystyle _{t_0}^t}𝑑t_1e^{\eta (tt_1)}\widehat{U}^{(0)}(t,t_1)`$ (15)
$`\times `$ $`[{\displaystyle \frac{d}{dt_1}}\sigma ^{(0)}(t_1)+{\displaystyle \frac{i}{\mathrm{}}}[\widehat{\overline{H}}^{(0)}(t_1),\sigma ^{(0)}(t_1)]`$
$`+`$ $`{\displaystyle \frac{i}{\mathrm{}}}[\widehat{\overline{H}}^{(1)}(t_1),\sigma (t_1)]]\widehat{U}^{(0)}^{}(t,t_1).`$
Since the Hamiltonian $`\widehat{\overline{H}}^{(0)}(t)`$ depends on time through the $`\widehat{Q}(t)`$, the time evolution operator $`\widehat{U}^{(0)}`$ is in general a time-ordered exponent

$$\widehat{U}^{(0)}(t,t_0)=T\mathrm{exp}\left[\frac{i}{\mathrm{}}_{t_0}^t𝑑t_1\widehat{\overline{H}}^{(0)}(t_1)\right].$$
(16)
To establish the time derivative $`d\sigma ^{(0)}(t_1)/dt_1`$ in Eq. (15), we recall from Eq. (9) that the relevant operator depends only implicitly on time through the averages $`\gamma _k(t_1)`$. Moreover by exploiting the transformation properties of a quantum Gaussian (see Appendix A), one finds that

$`{\displaystyle \frac{d}{dt}}\sigma ^{(0)}(t)`$ $`=`$ $`{\displaystyle \frac{d}{dt}}\gamma _k(t)_{\gamma _k}\sigma ^{(0)}(t)={\displaystyle \frac{i}{\mathrm{}}}[\widehat{\overline{H}}^{(0)}(t),\sigma ^{(0)}(t)]`$ (17)
$`+`$ $`{\displaystyle \frac{i}{\mathrm{}}}\mathrm{Tr}\left\{[\widehat{\overline{H}}^{(1)}(t),\widehat{\gamma }_k]\sigma (t)\right\}_{\gamma _k}\sigma ^{(0)}(t).`$
Using these relations we eliminate the time derivative of $`\sigma ^{(0)}(t_1)`$ in Eq. (15), to obtain the integral form of the statistical operator

$`\sigma (t)`$ $`=`$ $`\sigma ^{(0)}(t){\displaystyle \frac{i}{\mathrm{}}}{\displaystyle _{t_0}^t}𝑑t_1e^{\eta (tt_1)}\widehat{U}^{(0)}(t,t_1)`$ (18)
$`\times `$ $`[\text{Tr}\left\{[\widehat{\overline{H}}^{(1)}(t_1),\widehat{\gamma }_k]\sigma (t_1)\right\}_{\gamma _k}\sigma ^{(0)}(t_1)`$
$`+`$ $`[\widehat{\overline{H}}^{(1)}(t_1),\sigma (t_1)]]\widehat{U}^{(0)}^{}(t,t_1).`$
Since we are only interested in a weakly interacting gas, we seek a power series expansion in the interaction strength

$`\sigma (t)`$ $`=`$ $`\sigma ^{(0)}(t)+\sigma ^{(1)}(t)+\mathrm{}`$ (19)
$`\sigma ^{(1)}(t)`$ $`=`$ $`{\displaystyle \frac{i}{\mathrm{}}}{\displaystyle _{t_0}^t}𝑑t_1e^{\eta (tt_1)}\widehat{U}^{(0)}(t,t_1)`$ (20)
$`\times `$ $`[\text{Tr}\left\{[\widehat{\overline{H}}^{(1)}(t_1),\widehat{\gamma }_k]\sigma ^{(0)}(t_1)\right\}_{\gamma _k}\sigma ^{(0)}(t_1)`$
$`+`$ $`[\widehat{\overline{H}}^{(1)}(t_1),\sigma ^{(0)}(t_1)]]\widehat{U}^{(0)}^{}(t,t_1).`$
With this explicit expression for the statistical operator, the evaluation of the equation of motion (12) is straight forward and one obtains the quantum-Boltzmann equation

$`{\displaystyle \frac{d}{dt}}f(t)`$ $`=`$ $`[f]+[f]^{}.`$ (21)
Here, Wick’s theorem has been used to express the higher order averages in terms of the single particle ones. The kinetic operator $``$ consists of a reversible Hartree-Fock (HF) part $`_{\mathrm{HF}}`$ and a collisional quantum-Boltzmann contribution $`_{\mathrm{QB}}`$,

$`[f]`$ $`=`$ $`_{\text{HF}}[f]+_{\text{QB}}[f],`$
$`_{\text{HF}}[f]`$ $`=`$ $`{\displaystyle \frac{i}{\mathrm{}}}H_{\text{HF}}(t)f(t),`$
$`_{\text{QB}}[f]`$ $`=`$ $`\mathrm{\Gamma }_{ff(1+f)(1+f)}\mathrm{\Gamma }_{(1+f)(1+f)ff}.`$
$`H_{\mathrm{HF}}=\epsilon +2U_f`$ is the Hartree-Fock Hamiltonian, with $`U_f^{14}=2\varphi ^{1234}f_{32}`$ the self energy, while $`\mathrm{\Gamma }`$’s are the collision integrals given by

$`\mathrm{\Gamma }_{ABCD}^{15}`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{}}}{\displaystyle _{t_0}^t}dt_1[e^{\eta (tt_1)}\varphi ^{1234}\varphi ^{1^{\prime \prime }2^{\prime \prime }3^{\prime \prime }4^{\prime \prime }}`$ (22)
$`\times `$ $`𝒦_{1^{\prime \prime }1^{}}(t,t_1)𝒦_{2^{\prime \prime }2^{}}(t,t_1)𝒦_{3^{\prime \prime }3^{}}^{}(t,t_1)𝒦_{4^{\prime \prime }4^{}}^{}(t,t_1)`$
$`\times `$ $`A_{31^{}}(t_1)B_{42^{}}(t_1)C_{4^{}2}(t_1)D_{3^{}5}(t_1)].`$
The propagators $`𝒦(t,t_0)`$ are given by

$$𝒦(t,t_0)=T\mathrm{exp}\left[\frac{i}{\mathrm{}}_{t_0}^t𝑑t_1H_{\mathrm{HF}}(t_1)\right].$$
(23)
Therefore, unlike the collision terms of Ref. walser , the $`\mathrm{\Gamma }`$’s defined above depend on the past history of the system. Thus Eq. (21) represents a non-Markovian generalization of the kinetic equation previously derived in Ref. walser .

## III Conservation laws and Quasi-particle damping

The conserved quantities for a closed isolated system are the total energy $`E`$ and the total number $`N`$. These quantities then represent the constants of motion for the full kinetic equation. The total number operator can be represented as a linear combination of the relevant operators,

$$\widehat{N}=\widehat{a}_1^{}\widehat{a}_1^{},$$
(24)
and therefore the functional $`N(f)`$, representing the total number is given by

$$N(f)=\mathrm{Tr}\{f\}.$$
(25)
The kinetic equation for $`N(f)`$ can then be written as

$`{\displaystyle \frac{}{t}}N(f)`$ $`=`$ $`\mathrm{Tr}\left\{[f]+[f]^{}\right\}`$ (26)
$`=`$ $`\mathrm{Tr}\left\{_{HF}[f]+_{HF}[f]^{}\right\}`$
$`+`$ $`\mathrm{Tr}\left\{_{QB}[f]+_{QB}[f]^{}\right\}.`$
The first order term on the right-hand side involves the trace of a commutator and is trivially equal to zero. The $`\mathrm{\Gamma }`$’s associated with the second order terms have the following property due to the symmetries of $`\varphi `$,

$$\mathrm{Tr}\left\{\mathrm{\Gamma }_{ff(1+f)(1+f)}\right\}=\mathrm{Tr}\left\{\mathrm{\Gamma }_{(1+f)(1+f)ff}\right\}^{}.$$
(27)
As a result, the second order contribution in Eq. (26) can also be shown to be zero. Hence the total number $`N`$ is a constant of motion. The important point here is that the total number conservation is a result of the symmetries of $`\varphi `$ and does not dependent on the non-Markovian nature of the collision integral. Therefore a Markov approximation would leave this conservation law unchanged.

While the number conservation is a natural consequence of the self-consistency condition, the total energy conservation is not obvious as the Hamiltonian $`\widehat{H}`$ cannot be represented as a linear combination of the relevant operators. We start with writing the total energy functional $`E(f)`$ as a perturbative expansion in $`\varphi `$,

$`E(f)`$ $`=`$ $`\mathrm{Tr}\{\widehat{H}^{(0)}\sigma ^{(0)}\}`$ (28)
$`+`$ $`\left[\mathrm{Tr}\{\widehat{H}^{(1)}\sigma ^{(0)}\}+\mathrm{Tr}\{\widehat{H}^{(0)}\sigma ^{(1)}\}\right]`$
$`+`$ $`\left[\mathrm{Tr}\{\widehat{H}^{(1)}\sigma ^{(1)}\}+\mathrm{Tr}\{\widehat{H}^{(0)}\sigma ^{(2)}\}\right]+\mathrm{}`$
Using the self-consistency condition Eq. (11) for the relevant operators, one can show that the third and the fifth trace terms in the right-hand side of the above expression drop out. The kinetic equation for the energy functional $`E(f)`$ can then be written as

$`{\displaystyle \frac{}{t}}E(f)`$ $`=`$ $`\mathrm{Tr}\left\{\epsilon \dot{f}+U_f\dot{f}+U_{\dot{f}}f\right\}`$
$`+`$ $`{\displaystyle \frac{i}{2}}\mathrm{Tr}\left\{{\displaystyle \frac{}{t}}\left[\mathrm{\Gamma }_{ff(1+f)(1+f)}\mathrm{\Gamma }_{(1+f)(1+f)ff}\right]\right\}.`$
Again, we use the symmetry properties of $`\varphi `$, to write the simplified equation

$`{\displaystyle \frac{}{t}}E(f)`$ $`=`$ $`i{\displaystyle \frac{\eta }{2}}\mathrm{Tr}\left\{\mathrm{\Gamma }_{ff(1+f)(1+f)}\mathrm{\Gamma }_{(1+f)(1+f)ff}\right\}.`$
One can now see that the energy is conserved only in the $`\eta 0`$ limit. In this limit the kinetic equation (21) represents the Born approximation. A finite value of $`\eta `$ could then be thought of as resulting in additional terms that are beyond the Born approximation. Such terms model the duration of collision effects and quasi-particle damping. In principle if this effect is treated self-consistently, $`\eta `$ will be a time dependent function at least of order $`\varphi `$. This means that the rate of change of $`E`$ given by Eq. (LABEL:colle) is of the order $`\varphi ^3`$.

Also, if one is only interested in times grater than the correlation time $`\tau _{\mathrm{cor}}`$, the finiteness of $`\eta `$ allows us to extend the lower limit of the collision integral to $`\mathrm{}`$. Now we can approximate the $`f(t^{})`$ in the non-Markovian expression of $`\mathrm{\Gamma }`$ by its instantaneous value $`f(t)`$ to obtain the Markov form

$`\mathrm{\Gamma }_{ABCD}^{(m)^{15}}`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{}}}{\displaystyle _{\mathrm{}}^t}dt_1[e^{\eta (tt_1)}\varphi ^{1234}\varphi ^{1^{\prime \prime }2^{\prime \prime }3^{\prime \prime }4^{\prime \prime }}`$ (31)
$`\times `$ $`𝒦_{1^{\prime \prime }1^{}}(t,t_1)𝒦_{2^{\prime \prime }2^{}}(t,t_1)𝒦_{3^{\prime \prime }3^{}}^{}(t,t_1)𝒦_{4^{\prime \prime }4^{}}^{}(t,t_1)`$
$`\times `$ $`A_{31^{}}(t)B_{42^{}}(t)C_{4^{}2}(t)D_{3^{}5}(t)]`$
One can verify that the above Markov form results in an energy conservation, up to the most significant order given by

$`{\displaystyle \frac{}{t}}E(f)`$ $`=`$ $`\mathrm{Tr}\left\{H_{\mathrm{HF}}\left(_{\mathrm{QB}}^{(m)}[f]+_{\mathrm{QB}}^{(m)}[f]^{}\right)\right\}.`$ (32)
If we compare this with the expression for the correlation energy $`E_{\mathrm{cor}}`$ in Refs. morozov ; morozov2 , the right hand side of Eq. (32) is exactly $`E_{\mathrm{cor}}/t`$. This is not surprising and can be understood more intuitively by writing the collision integral as a sum of two contributions: correlation and collision

$$\mathrm{\Gamma }=\mathrm{\Gamma }_{\mathrm{cor}}+\mathrm{\Gamma }_{\mathrm{col}}=_{\mathrm{}}^0\mathrm{}𝑑t^{}+_0^t\mathrm{}𝑑t^{}.$$
(33)
For a finite $`\eta `$, the $`\mathrm{\Gamma }_{\mathrm{cor}}`$ contribution decays to zero as $`e^{\eta t}`$. Therefore the decaying correlation energy $`E_{\mathrm{cor}}`$ associated with this part shows up in the rate of change of the total energy $`E`$. In the Born approximation $`\mathrm{\Gamma }_{\mathrm{cor}}`$ is constant because $`\eta 0`$.

The exponential damping results in a widened delta function and therefore the rate of change of $`E(f)`$ Eq. (32) can be shown to be of order $`\eta \mathrm{\Gamma }^{(m)}`$. Thus, by including terms beyond the Born approximation, we have obtained a collision integral that is Markovian and still conserves energy up to $`\varphi ^2`$ order. Now if we assume that the equilibration time is of the order $`1/\mathrm{\Gamma }^{(m)}`$ , the total change in energy $`\mathrm{\Delta }E`$ is therefore

$`\mathrm{\Delta }E`$ $`=`$ $`E(f^{\mathrm{eq}})E(f^{\mathrm{in}})`$ (34)
$``$ $`\eta \mathrm{\Gamma }^{(m)}\times \left({\displaystyle \frac{1}{\mathrm{\Gamma }^{(m)}}}\right)=\eta ,`$
where $`f^{\mathrm{in}}`$ and $`f^{\mathrm{eq}}`$ are the initial and equilibrium distributions, respectively.

The importance of the damping term becomes more obvious when one attempts to solve the kinetic equation numerically. One no longer has to worry about the $`\eta 0`$ limit in the Born-Markov approximation, which for a finite system with discrete levels will result only in exchange collisions and hence no equilibration.

## IV  Application to a dilute Bose gas in a spherical box trap

In the previous section, we introduced the general methods and concepts to describe a weakly interacting Bose gas under nonequilibrium conditions. We will now apply these to a simple model of a typical <sup>87</sup>Rb experiment, as realized by many laboratories around the world, for example Refs. cornell ; jin96 ; jin97 . The physical parameters are usually quoted in the natural units for a harmonic oscillator trapping potential, i. e., the angular frequency $`\omega =2\pi \mathrm{\hspace{0.17em}200}\text{Hz}`$, the atomic mass $`m_{87}=86.9092\text{amu}`$, the ground state size $`a_\text{H}=[\mathrm{}/(\omega m_{87})]^{1/2}=763\text{nm}`$, and the s-wave scattering length $`a_\text{S}=5.82\text{nm}`$.

However, in the present article, we do not pursue the usual harmonic confinement, but rather explore the properties of a radial box as a particle trap. This choice is motivated by previous studies of the self-consistent Hartree-Fock single-particle states walser2 . As soon as repulsive mean-field potentials are added to the bare harmonic trapping potentials, the corresponding eigenstates widen in size and look remarkably close to the eigenstates of a box, provided the spatial extensions of the box is chosen appropriately. In particular, we pick a box of radius $`R=1000a_\text{S}=5.82\mu \text{m}`$.

Thus, our model is represented by a spherical trap with box potential given by

$$V_{\mathrm{ext}}(r)=\{\begin{array}{cc}0,& r<R\hfill \\ \mathrm{},& rR,\hfill \end{array}$$
(35)
The eigenfunctions are a product of spherical Bessel functions $`j_{(l)}`$ and spherical harmonics $`Y_{(lm)}`$:

$$\psi _{(nlm)}(r,\theta ,\varphi )=\{\begin{array}{cc}𝒩_{(nl)}j_{(l)}(rk_{(nl)})Y_{(lm)}(\theta ,\varphi ),\hfill & r<R\hfill \\ 0,\hfill & rR.\hfill \end{array}$$
(36)
A normalization constant is represented by $`𝒩_{(nl)}`$. The eigenenergies are given in terms of the wave-vectors $`k_{(nl)}`$, which can be obtained from the $`n`$-th nodes of the spherical Bessel-functions of angular momentum $`l`$,

$`\epsilon _{(nl)}`$ $`=`$ $`(Rk_{(nl)}/\pi )^2\epsilon _0.`$ (37)
Here, $`\epsilon _0`$ represents the ground state energy and defines the energy scale of the problem

$`\epsilon _0`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2\pi ^2}{2mR^2}}.`$ (38)
All the physical parameters are scaled with respect to this energy unit $`\epsilon _0`$ and the radius of the box $`R`$.

For simplicity, we assume all the atoms to be in the $`l=0`$ state initially. One therefore needs to consider only the $`l=0`$ manifold and get for the normalization constant $`𝒩_{(n0)}=n(\pi /(2R^3))^{1/2}`$ and radial wave-vector $`Rk_{(n0)}=n\pi `$. If the cloud is relatively cold, then most of the population resides in the lowest few energy states. Therefore, we can also limit the number of radial modes $`1nn_{\text{max}}`$. For the present case, we take $`n_{\text{max}}=5`$. Obviously, all these simplification reduce the number of degrees of freedom significantly and thus we are able to study certain aspects of the nonequilibrium dynamics of the trapped Bose gas in great detail.

With the above definitions, the bare single particle box Hamiltonian is given by

$$\epsilon ^{pq}=\delta _{pq}q^2(\text{no implicit sum over}q).$$
(39)
The interaction part $`\widehat{H}^{(1)}`$ involves the matrix elements of the interaction potential defined by Eq. (6)

$`\varphi ^{pqrs}`$ $`=`$ $`{\displaystyle \frac{4a_\text{S}}{\pi }}{\displaystyle _0^\pi }\mathrm{sin}(px)\mathrm{sin}(qx)\mathrm{sin}(rx)\mathrm{sin}(sx){\displaystyle \frac{dx}{x^2}},`$
which in general have to be computed numerically. Interestingly, in the case of a spherical box this integral can be evaluated analytically and simplified to a finite sum of sine integrals and cosine functions (see Appendix B for details).

The nonequilibrium state of the above system is represented by the single particle distribution function $`f`$ with time dependence given by Eq. (21). Neglecting the second order collision terms, the first order evolution is governed by the Hartree-Fock Hamiltonian $`H_{\mathrm{HF}}`$ and is given by

$$\frac{}{t}f=_{HF}[f]+_{HF}[f]^{}.$$
(41)
The energy eigenstates of the interacting system are therefore shifted from the bare box states due to the self energy effect. These shifts can be significant depending on the total particle number. This is clear from Fig. 1 where we plot the eigenenergies as a function of $`f_{11}`$, the total particle number in the box ground state (with all other $`f_{ij}`$’s equal to zero).

Now note that the time dependent Hartree-Fock equation (41) for the density matrix, $`f`$, is nonlinear and hence we seek a self consistent solution such that

$$f=\underset{i}{}P(\epsilon _i)|\epsilon _i\epsilon _i|,$$
(42)
where $`H_{\mathrm{HF}}|\epsilon _i=\epsilon _i|\epsilon _i`$, and for a Bose-Einstein distribution $`P(\epsilon )`$ is given by

$$P(\epsilon )=\frac{1}{\mathrm{exp}((\epsilon \mu )/k_BT)1}.$$
(43)
For a given total particle number, $`N`$, and temperature, $`\beta =1/k_BT`$, a self-consistent chemical potential, $`\mu (\beta ,N)`$, and hence a self-consistent Bose-Einstein distribution is obtained. For example, let us consider three different total particle numbers $`N=10,100,500`$ at two different temperatures, $`\beta =1/k_BT=0.01,0.5`$, corresponding to hot and cold clouds, respectively. In Figs. 2 and 3, we plot the self-consistent solution and the self-energy density in the position space representation. We see that the self-energy density is proportional to the number density only near the center of the trap and drops off faster with increasing radius. This can be attributed to the restricted number of basis states used in our calculation and effectively gives a finite range to the two body potential $`V_{\mathrm{bin}}`$.

Up to first order, the $`f`$ equation is totally reversible. The inclusion of the second order terms (collisions) break the reversibility, and therefore represents a relaxation of the system from some initial state to a final equilibrium state. Here we will be using the Markov form (31) for the collision integral. We take the Hartree-Fock self-consistent state for the initial condition.

In the previous section, we interpreted the function $`\mathrm{exp}(\eta \tau )`$ to account for duration of collision effects and quasi-particle damping. But the exponential form was originally introduced to break the time reversal symmetry and it has the correct long time behavior. However, an exponential damping will result in a Lorenzian line shape for the final equilibrium distribution. Due to the long-reaching wings of the Lorenzian curve, in the Markov limit, off-the-energy shell collisions get weighted strongly. To seek an improved damping function that will have the correct short and long time behavior, we use the equivalence of kinetic theories based on the Green’s function approach kadanoff  and the non-equilibrium statistical operator method zubarev  as shown in Ref. wachter . The behavior of the retarded Green’s function $`g(t,t_1)`$ for very large and very small time scales is given by

$$g(t,t_1)\{\begin{array}{cc}e^{\eta (tt_1)}\hfill & tt_1\tau _{\mathrm{cor}}\hfill \\ e^{\eta (tt_1)^2}\hfill & tt_1\tau _{\mathrm{cor}}.\hfill \end{array}$$
(44)
Therefore, behavior over the intermediate time scale will be best represented by an interpolating function. This is also true for the damping function. From Fig. 4 we see that the function

$$(t,t_1)=1/\mathrm{cosh}(\eta (tt_1)),$$
(45)
has exactly this behavior and therefore represents a better choice than the exponential form.

With either choice of the damping function, for a particular value of the parameter $`\eta `$, a time propagation results in a self-consistent steady-state solution, $`f^{\mathrm{eq}}`$. Figure 5 shows such a time evolution for $`N=500`$ particles with an initial temperature corresponding to $`\beta =0.01`$. We have chose $`\eta `$ to be of the order $`\mathrm{Re}[\mathrm{\Gamma }]2.3`$. As mentioned earlier we see that the exponential damping function results in significant transfer of population to the excited states. This effect is less with the hyperbolic secant damping function. We get different steady-state solutions because different damping functions correspond to different initial correlations. Also a plot of $`\theta (f_{nn})\mathrm{ln}(1/f_{nn}+1)`$ vs the Hartree-Fock energies $`\epsilon _n^{HF}`$, shown in Fig. 6 shows that the $`f^{\mathrm{eq}}`$ is very close to a Bose-Einstein distribution. The slope, which represents the self-consistent value of $`\beta `$, shows that the change in temperature is far greater for the exponential damping function compared to the hyperbolic secant case, because off-the-energy shell effects are larger as explained previously.

Thus we have obtained a self-consistent steady-state solution to the second order kinetic equation. We emphasize here that the steady-state solution is a result of the real-time non-equilibrium evolution of the system. Deriving such an equilibrium solution, whose absolute value is plotted in Fig. 7, is a prerequisite step to the study of collective modes and damping rates of a dilute gas.

## V Real time response to perturbation

The properties of the equilibrium solution $`f^{\text{eq}}`$ exhibit the expected characteristics of a Bose-Einstein distribution. In order to verify the stability of this solution and to study the damping rates of the collective excitations, we will now examine the real-time response of the system to a perturbation. First, we will outline the linear response theory and discuss the structure of the modes, their frequencies, and the life-time of the excitations. Subsequently, we will use these modes to initially prepare the system and to evolve the full nonlinear quantum kinetic equation towards equilibrium.

One of the fundamental properties of the quantum kinetic equation (21)

$`{\displaystyle \frac{d}{dt}}f(t)`$ $`=`$ $`[f]+[f]^{},`$ (46)
is its Hermitian structure. Thus, if we prepare a physical state initially, it will remain Hermitian with $`f(t)=f(t)^{}`$, indefinitely. We will now consider a weak perturbation of an equilibrium state,

$`f(t)`$ $`=`$ $`f^{\text{eq}}+\delta f(t),`$ (47)
and calculate the first order response of the system. As usual, we want to assume that we can decompose a general perturbation into fundamental damped and/or oscillatory eigenmodes of the system. Therefore, such a specific perturbation can be parameterized as

$`\delta f(t)`$ $`=`$ $`e^{i\omega t}\delta f_\omega ^{(+)}+\text{h.c.},`$ (48)
$`\delta f_\omega ^{(+)}`$ $`=`$ $`\delta f_\omega ^{(c)}+i\delta f_\omega ^{(s)},`$ (49)
where $`\delta f_\omega ^{(+)}`$ denotes a positive frequency amplitude. It turns out to be useful to decompose it into quadrature components

$`\delta f_\omega ^{(c)}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left[\delta f_{\omega }^{(+)}{}_{}{}^{}+\delta f_\omega ^{(+)}\right],`$ (50)
$`\delta f_\omega ^{(s)}`$ $`=`$ $`{\displaystyle \frac{1}{2i}}\left[\delta f_{\omega }^{(+)}{}_{}{}^{}\delta f_\omega ^{(+)}\right],`$ (51)
and evaluate the kinetic operator $`[f]`$ only for such Hermitian arguments. From a Taylor series expansion of the kinetic operator around the equilibrium distribution, one obtains finally the linear response equations for the fundamental modes

$`(i\omega )\delta f_\omega ^{(+)}`$ $`=`$ $`^{(1)}[\delta f_\omega ^{(c)}]+i^{(1)}[\delta f_\omega ^{(s)}],`$ (52)
$`+`$ $`^{(1)}[\delta f_\omega ^{(c)}]^{}+i^{(1)}[\delta f_\omega ^{(s)}]^{},`$
where we have defined a linear response operator through an appropriate centered difference limit

$`^{(1)}[\delta f]`$ $`=`$ $`\underset{\lambda 0}{lim}{\displaystyle \frac{[f^{\text{eq}}+\lambda \delta f][f^{\text{eq}}\lambda \delta f]}{2\lambda }}.`$ (53)
We solve Eq. (52) as an eigenvalue problem. In general, the eigenvalues are complex, with frequency and damping rate given by the real and the imaginary parts respectively. The eigenvalues appear as complex conjugate pairs. The eigenmodes corresponding to non-zero eigenvalues are Hermitian conjugates of each other and are traceless with normalization given by

$$\mathrm{Tr}\{\delta f_\omega ^{(+)}\delta f_\omega ^{(+)}{}_{}{}^{}\}=1.$$
(54)
The physical linear response mode is given by the quadrature components $`\delta f_\omega ^s`$ and $`\delta f_\omega ^c`$. There also exists a zero mode that has nonvanishing trace. The damping rates are all negative, thus confirming the stability of the collective modes. In Fig. 8, we plot the positive frequency eigenvalues. The dotted lines correspond to the difference frequencies of the Hamiltonian $`H_{HF}`$.

It is interesting to see how these different modes evolve in real time. For this we use the equilibrium distribution obtained in the previous section and perturb it with one of the quadrature components

$`ff^{\mathrm{eq}}+\lambda \delta f_\omega ^{(s)},`$ (55)
where $`\lambda =0.2`$ determines the smallness of the perturbation. In particular we will consider the modes labeled by (a), (b), (c) and (d) in Fig. 8.

The real-time response is shown by plotting the off-diagonal matrix element $`f_{12}`$ of the single particle density matrix in the box basis as shown in Figs. 9 and 10. In Fig. 11, we plot the change in the total energy $`\mathrm{\Delta }E=E(f)E(f^{\mathrm{eq}})`$ as a function of time. In cases (a) and (b) we see that the $`\mathrm{\Delta }E`$ oscillates about zero and eventually goes to zero. This is expected because such a perturbation tends to create coherences, resulting in an energy change by the amount of $`E_{\mathrm{coh}}`$ (coherence energy) that would eventually decay down to zero. Similar damped behavior is

observed for case (d). Such an oscillatory damped behavior of the total energy could be attributed to the Markov approximation in the collision integral Eq. (31). On the other hand, perturbations of the kind (c) increase the total particle number by the amount $`\delta N=\mathrm{Tr}\{\lambda \delta f_\omega ^{(s)}\}`$ and hence result in a finite change in total energy.

## VI Conclusion

A non-Markovian version of the quantum kinetic theory is derived using the prescription of a nonequilibrium statistical operator method as outlined in Ref. zubarev . This theory is shown to conserve energy in the $`\eta 0`$ limit. Inclusion of quasi-particle damping results in a description beyond the Born approximation with energy conservation to $`\varphi ^2`$ order even in the Markov limit. To obtain collision integrals that involve quasi-particle damping and conserves energy precisely, one will have to calculate the $`T`$-matrix in the full collision operator keeping terms of all orders in the interaction.

We applied the generalized second order kinetic theory to the nonhomogeneous dilute Bose gas confined in a spherical box to numerically study the full non-equilibrium evolution of the system towards equilibrium. The self-consistent distribution $`f^{\mathrm{eq}}`$ thus obtained is very close to the Bose-Einstein distribution as shown in Fig. 6. We also observe a significant Hartree-Fock self energy shift which depends on the single particle distribution function $`f`$. The form of the damping function is important in determining the line shape. Particularly, the function with a $`1/\mathrm{cosh}`$ type of behavior is found to be appropriate and gives improved energy conservation due to smaller initial correlation effect.

The importance of such a real-time calculation is apparent from the full real-time response calculation, where we have calculated the damping rates and frequencies. These damping rates correspond to a shorter time scale compared to the equilibration time scale, which depend on rates in and out of the various levels.

This simple model of a spherical trap can be easily extended to a more realistic situation of a harmonic trap. As in Refs. walser ; kane ; martin , the condensed component can be easily included by introducing a symmetry broken mean-field, $`\alpha _i=\widehat{a}_i`$, as one of the relevant observables and Hartree-Fock-Bogoliubov quasi-particle excitations. Such a calculation will allow us to make experimentally verifiable predictions of damping rates of collective excitations. One can also explore the possibility of including a time dependent potential or an external force term to selectively excite one or more of the collective modes. These possibilities will be explored in a forthcoming article wachter2 .

###### Acknowledgements.

We are extremely thankful to J. Cooper for the numerous and invaluable discussions. S.B. and M.H. acknowledge the financial support from the U. S. Department of Energy, Office of Basic Energy Sciences via the Chemical Sciences, Geosciences and Biosciences Division. R.W. gratefully receives financial support from the Austrian Academy of Sciences through an APART grant.

## Appendix A Reference distribution

The reference distribution $`\sigma _{\{f\}}^{(0)}`$ of Eq. (9) is parameterized through its expectation values in Eq. (11). From the structure of this quantum Gaussian operator, it follows that

$$\sigma _{\{𝒦(t,t_0)f𝒦^{}(t,t_0)\}}^{(0)}=\widehat{U}^{}(t,t_0)\sigma _{\{f\}}^{(0)}\widehat{U}(t,t_0).$$
(56)
In the above, $`\widehat{U}`$ represents the single particle propagator of Eq. (16) acting in many-particle Fock-space and $`𝒦`$ is the corresponding single particle Hilbert-space propagator of Eq. (23). This condition implies that

$$_{\gamma _k}\sigma ^{(0)}\text{Tr}\{[\widehat{H}^{(0)},\widehat{\gamma }_k]\}=[\widehat{H}^{(0)},\sigma ^{(0)}]$$
(57)
and was used to obtain Eq. (17).

## Appendix B Matrix elements

The matrix element of Eq. (IV) can be evaluated easily by expanding the $`\mathrm{sin}`$ functions into co- and counter-propagating complex exponents and by an additional partial integration. This results in eight separate terms, i.e.,

$`\varphi ^{ijkl}/a_\text{S}`$ $`=`$ $`{\displaystyle \frac{4}{\pi }}{\displaystyle _0^\pi }\mathrm{sin}(ix)\mathrm{sin}(jx)\mathrm{sin}(kx)\mathrm{sin}(lx){\displaystyle \frac{dx}{x^2}},`$ (58)
$`=`$ $`F(i+j+kl)+F(i+jk+l)`$
$`+`$ $`F(ij+k+l)+F(ijkl)`$
$``$ $`F(i+jkl)F(ij+kl)`$
$``$ $`F(ijk+l)F(i+j+k+l),`$
where

$`F(n)`$ $`=`$ $`{\displaystyle \frac{1}{2\pi ^2}}[\mathrm{cos}(n\pi )+n\pi \text{Si}(n\pi )],`$ (59)
$`\mathrm{Si}(z)`$ $`=`$ $`{\displaystyle _0^z}{\displaystyle \frac{\mathrm{sin}(t)}{t}}𝑑t.`$ (60)
An asymptotic expansion of the sine integral leads to the following approximation that is correct at the $`1\%`$ level, i.e.,

$`F(0)`$ $`=`$ $`{\displaystyle \frac{1}{2\pi ^2}},`$ (61)
$`F(n>0)`$ $``$ $`{\displaystyle \frac{1}{2\pi ^2}}[{\displaystyle \frac{\pi }{2}}|n\pi |{\displaystyle \frac{\mathrm{sin}(n\pi )}{n\pi }}+{\displaystyle \frac{2\mathrm{cos}(n\pi )}{(n\pi )^2}}].`$