text:
# Floquet scattering theory of quantum pumps

## I Introduction

Quantum charge pumping  is presently of considerable interest. An experiment by Switkes et al.  demonstrated that a phase coherent mesoscopic system subjected to a cyclic two parameter perturbation can produce a directed current.

Coherent quantum pumping is a consequence of the interference of energetically different traversal paths made possible by an oscillating scatterer. The ratio of the oscillation frequency $`\omega `$ to the inverse time taken for carriers to traverse the sample $`\tau _T^1`$ defines the operational regime of a pump . Brouwer  gave an elegant formulation of adiabatic ($`\omega \tau _T^1`$) quantum pumping that is based on the scattering matrix approach to low frequency ac transport in phase coherent mesoscopic systems developed by Büttiker, Thomas, and Prêtre . This approach leads naturally to a geometrical description of adiabatic quantum pumping . The theory predicts that the charge pumped during a cycle depends on the area enclosed by the path in the scattering matrix parameter space. A less formal but more physical picture of an adiabatic quantum pump appeals to both quantum mechanical interference and photon assisted transport. . The same processes are important for a nonadiabatic ($`\omega \tau _T^1`$) pump . These discussions emphasize the energetics of the carrier traversal process.

It is the purpose of this work to develop a theory that permits the description of both adiabatic and nonadiabatic regimes on the same footing and allows a simple physical interpretation. To this end we extend the approach of Ref. to the case of large frequencies and large pumping amplitudes. We apply the Floquet scattering theory  which deals with the scattering matrix dependent on two energies (incident and outgoing). This approach leads to expressions for the quantities of interest in terms of the side bands of particles exiting the pump. The side bands  correspond to particles which have gained or lost one or several modulation quanta $`\mathrm{}\omega `$. This approach is complementary to discussions based on the scattering matrix dependent on two times

The paper is organized as follows. In Sec.II the general approach to the kinetics of quantum pumps based on the Floquet scattering theory is presented. In Sec.III we apply the general results to the adiabatic case. In Sec.IV we calculate the Floquet scattering matrix for a particular model - an oscillating double barrier potential - and present the results of numerical calculations of the pumped charge and the heat currents in both adiabatic and nonadiabatic regimes. We conclude in Sec.V.

## II General approach

We consider scattering  of an incoming flow of electrons with energy $`E`$ at a scatterer that oscillates in time with frequency $`\omega `$.

During the interaction with the oscillating scatterer  electrons can gain or lose energy quanta $`\mathrm{}\omega `$. Hence the outgoing state is characterized by the set of energies $`E_n`$, $`n=0,\pm 1,\pm 2,\mathrm{}`$,

$$E_n=E+n\mathrm{}\omega .$$
(1)
This is a Floquet state.

According to the Floquet theorem the energy ladder Eq.(1) gives the full set of possible energies for outgoing particles (see e.g., ). Thus to describe scattering due to an oscillating scatterer we can use the Floquet scattering matrix $`\widehat{S}_F`$. The matrix element $`S_{F,\alpha \beta }(E_n,E)`$ is the quantum mechanical amplitude for an electron with energy $`E`$ entering the scatterer through lead $`\beta `$ to leave the scatterer through lead $`\alpha `$ having absorbed ($`n>0`$) or emitted ($`n<0`$) energy quanta $`|n|\mathrm{}\omega `$. The Greek letters $`\alpha `$, $`\beta `$, number the leads connecting the sample to $`N_r`$ reservoirs.

We remark that the negative values $`E_n<0`$ correspond to bound states near the oscillating scatterer. These states influence scattering into the propagating ($`E_n>0`$) states but they do not directly contribute to the current.

Current conservation implies that the submatrix $`\widehat{S}_F^{(p)}`$ of the Floquet scattering matrix (corresponding to propagating modes only) is a unitary matrix

$$\widehat{S}_F^{(p)}\widehat{S}_F^{(p)}=\widehat{S}_F^{(p)}\widehat{S}_F^{(p)}=\widehat{I}.$$
(2)
In particular, if a current with flux $`1`$ and energy $`E`$ enters the scatterer through lead $`\beta `$ then current conservation implies

$$\underset{\alpha }{}\underset{E_n>0}{}|S_{F,\alpha \beta }(E_n,E)|^2=1.$$
(3)
Another useful condition follows from the fact that if all incoming propagating ($`E_n>0`$) channels are full then each outgoing channel has also to be full:

$$\underset{\beta }{}\underset{E_n>0}{}|S_{F,\alpha \beta }(E,E_n)|^2=1.$$
(4)
Note that usually the Floquet energy $`E`$ is determined within the interval $`0E<\mathrm{}\omega `$. However, for our problem, it is convenient not to reduce the discrete set of $`E_n`$ to this interval and to keep $`E`$ as the actual energy of incident (or outgoing) particles.

Because of Eq.(2) we can express the annihilation operator $`\widehat{b}`$ for outgoing particles in the lead $`\alpha `$ in terms of annihilation operators $`\widehat{a}`$ for incoming particles in leads $`\beta =1,2,\mathrm{},N_r`$ as follows

$$\widehat{b}_\alpha (E)=\underset{\beta }{}\underset{E_n>0}{}S_{F,\alpha \beta }(E,E_n)\widehat{a}_\beta (E_n).$$
(5)
The operators $`\widehat{a}_\alpha (E)`$ obey the following anticommutation relations

$$[\widehat{a}_\alpha ^{}(E),\widehat{a}_\beta (E^{})]=\delta _{\alpha \beta }\delta (EE^{}).$$

Using Eqs.(2) and (5) we see that the operators $`\widehat{b}_\alpha (E)`$ obey the same relations.

Note that above expressions correspond to single (transverse) channel leads and spinless electrons. For the case of many-channel leads each lead index ($`\alpha `$, $`\beta `$, etc.) includes a transverse channel index and any repeating lead index implies implicitly a summation over all the transverse channels in the lead. Similarly an electron spin can be taken into account.

Now we calculate the distribution function $`f^{(out)}(E)_\alpha `$ $`=`$ $`<\widehat{b}_\alpha ^{}(E)\widehat{b}_\alpha (E)>`$ for electrons leaving the scatterer through the lead $`\alpha `$. Here $`<\mathrm{}>`$ means quantum-statistical averaging. Taking into account Eq.(5) we obtain

$$\begin{array}{c}f_\alpha ^{(out)}(E)=\underset{\beta }{}\underset{E_n>0}{}|S_{F,\alpha \beta }(E,E_n)|^2\\ \\ \times f_\beta ^{(in)}(E_n).\end{array}$$
(6)
Here $`f_\beta ^{(in)}(E_n)`$ $`=`$ $`<\widehat{a}_\beta ^{}(E)\widehat{a}_\beta (E)>`$ is the distribution function for electrons entering the scatterer through lead $`\beta `$.

### A Directed charge currents

Using the distribution function $`f_\alpha ^{(out)}(E)`$ for outgoing particles and $`f_\alpha ^{(in)}(E)`$ for incoming ones we can find the directed current $`I_\alpha `$ in the lead $`\alpha `$ far from the scatterer

$$I_\alpha =\frac{e}{h}\underset{0}{\overset{\mathrm{}}{}}𝑑E\left\{f_\alpha ^{(out)}(E)f_\alpha ^{(in)}(E)\right\}.$$
(7)
The current directed from the scatterer towards the reservoir is positive by definition. Substituting Eq.(6) into the above equation, using Eq.(4), and making the shift $`EEn\mathrm{}\omega `$ we find

$$\begin{array}{c}\hfill I_\alpha =\frac{e}{h}\underset{0}{\overset{\mathrm{}}{}}𝑑E\underset{\beta }{}\underset{E_n>0}{}|S_{F,\alpha \beta }(E_n,E)|^2\\ \\ \hfill \times \left(f_\beta ^{(in)}(E)f_\alpha ^{(in)}(E_n)\right).\end{array}$$
(8)
Here $`_{E_n>0}`$ means a sum over those $`n`$ (positive and negative) for which $`E_n=E+n\mathrm{}\omega >0`$.

Another useful representation for the directed current can be obtained if we use Eq.(3) and make the shift $`EEn\mathrm{}\omega `$ in $`f_\alpha ^{(out)}(E)`$ in Eq.(7). As a result we obtain

$$\begin{array}{c}I_\alpha =\frac{e}{h}\underset{0}{\overset{\mathrm{}}{}}𝑑E\underset{\beta \alpha }{}\underset{E_n>0}{}\hfill \\ \\ \{|S_{F,\alpha \beta }(E_n,E)|^2f_\beta ^{(in)}(E)\hfill \\ \\ |S_{F,\beta \alpha }(E_n,E)|^2f_\alpha ^{(in)}(E)\}.\hfill \end{array}$$
(9)
¿From this expression for the directed current it follows that only transmission $`\alpha \beta `$ (not reflection $`\alpha =\beta `$) contributes to the current. In addition Eq.(9) can help us to consider the effect of time reversal symmetry (TRS) on the pumped current.

On the one hand, the time reversal $`tt`$ (TR) interchanges incoming and outgoing channels

$$\left[S_{F,\alpha \beta }(E_n,E)\right]^{(TR)}=S_{F,\beta \alpha }(E,E_n).$$

Hence if the TRS is present then Eq.(9) reads

$$\begin{array}{c}I_\alpha ^{(TRS)}=\frac{e}{h}\underset{0}{\overset{\mathrm{}}{}}𝑑Ef_0(E)\underset{\beta \alpha }{}\underset{E_n>0}{}\hfill \\ \\ \left(|S_{F,\alpha \beta }(E_n,E)|^2|S_{F,\alpha \beta }(E,E_n)|^2\right).\hfill \end{array}$$
(10)
In the above equation (in accordance with the usual pump setup) we suppose that incoming electrons in all the channels are described by the same Fermi distribution function with temperature $`T`$ and electrochemical potential $`\mu `$:

$$f_\alpha ^{(in)}(E)=f_0(E)\frac{1}{1+\mathrm{exp}\left(\frac{E\mu }{k_BT}\right)}.$$

Generally, the Floquet scattering matrix elements for transmission with incident energy $`E`$ to $`E_n`$ is not equal to the transmission from $`E_n`$ to $`E`$,

$$S_{F,\alpha \beta }(E_n,E)S_{F,\alpha \beta }(E,E_n).$$

¿From this we can conclude that even a pump with TRS can generate a directed current. If these two scattering amplitudes are not equal, there exists the possibility of empty states deep below the Fermi surface . In this case interaction and inelastic effects  can be expected to be especially important.

However, if the scattering matrix is energy independent on the scale of the order of $`\mathrm{}\omega `$

$$S_{F,\alpha \beta }(E_n,E)S_{F,\alpha \beta }(E),$$

then the scatterer with TRS can not produce a dc current. The last circumstance is especially important for the adiabatic case $`\omega 0`$, since the adiabatic scattering matrix Eq.(16) satisfies the above condition.

On the other hand, in general, the nonstationary case is without TRS and the oscillating scatterer can, in principle, generate a current. In particular, if two parameters of the scatterer oscillate in time with the same frequency $`\omega `$ but with phase lag $`\mathrm{\Delta }\phi `$ then the time reversal implies the reversal of the sign of $`\mathrm{\Delta }\phi `$. Thus for a pump with $`\mathrm{\Delta }\phi 0`$ the scattering problem is without TRS (the time reversal symmetry is dynamically broken) and such a scatterer can produce a dc current. This effect was used for generating a dc current in adiabatic quantum pumps .

### B Directed heat currents

By analogy with Eq.(7) we find the directed heat current $`I_{E,\alpha }`$ flowing in lead $`\alpha `$ away from the scatterer  (we suppose that all the reservoirs are at the same electrochemical potential $`\mu `$)

$$\begin{array}{c}I_{E,\alpha }=\frac{1}{h}\underset{0}{\overset{\mathrm{}}{}}𝑑E(E\mu )\hfill \\ \\ \times \left\{f_\alpha ^{(out)}(E)f_\alpha ^{(in)}(E)\right\}.\hfill \end{array}$$
(11)
With the distribution function for outgoing particles $`f^{(out)}`$ from Eq.(6) we get

$$\begin{array}{c}I_{E,\alpha }=\frac{1}{h}\underset{0}{\overset{\mathrm{}}{}}𝑑E\underset{\beta }{}\underset{E_n>0}{}(E_n\mu )\hfill \\ \\ \times |S_{F,\alpha \beta }(E_n,E)|^2\left(f_\beta ^{(in)}(E)f_\alpha ^{(in)}(E_n)\right).\hfill \end{array}$$
(12)
Note that if all the reservoirs are at the same macroscopic conditions (electrochemical potential, temperature, etc.) then the heat flow $`I_{E,\alpha }`$ (at any lead $`\alpha =1,2,\mathrm{},N_r`$) is directed from the scatterer to the reservoir . That differs strongly from the charge current $`I_\alpha `$ given by Eq.(8) which, if it exists, can be directed either from the reservoir to the scatterer (at some lead) or vice versa (at another lead).

## III Adiabatic approximation

In this section we use the above formalism to investigate the limit of adiabatic scattering.

The general physical notion of adiabaticity applied to the scattering problem of interest here is as follows: Let us suppose that the time independent problem is described by the scattering matrix $`\widehat{S}_0(E,X_1,\mathrm{},X_N)`$ which depends on the energy $`E`$ of incident electrons and a set of parameters $`X_i`$, $`i=1,2,\mathrm{},N`$. Next assume that the parameters $`X_i`$ vary in time: $`X_i=X_i(t)`$. Then we can say that the nonstationary scattering problem is adiabatic if scattering of particles incident with energy $`E`$ can be described via the scattering matrix $`\widehat{S}_0`$ with time-dependent parameters $`X_i(t)`$

$$\widehat{S}_{ad}(E,t)=\widehat{S}_0(E,X_1(t),\mathrm{},X_N(t)).$$
(13)
This approximation is adequate if the scattering matrix changes only a little while an electron interacts with the scatterer. In other words, the characteristic time scale for the change of parameters $`X_i`$ is much larger than the traversal time $`\tau _T`$.

An analogous criterion can be formulated concerning the energy dependence of the scattering matrix. To this end we consider the adiabatic problem when the parameters change periodically in time: $`X_i(t)=X_i(t+2\pi /\omega )`$. In this case we can expand the adiabatic scattering matrix into the Fourier series

$$\widehat{S}_{ad}(E,t)=\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}\widehat{S}_{0,n}(E)e^{in\omega t},$$
(14)
where

$$\begin{array}{c}\widehat{S}_{0,n}(E)=\frac{\omega }{2\pi }\underset{0}{\overset{2\pi /\omega }{}}𝑑te^{in\omega t}\hfill \\ \\ \times \widehat{S}_0(E,X_1(t),\mathrm{},X_N(t)).\hfill \end{array}$$
(15)
The Fourier harmonics $`S_{0,n}`$ define the amplitudes of side bands for particles traversing the adiabatically oscillating scatterer with initial incident energy $`E`$. Thus we can construct the adiabatic Floquet scattering matrix as follows

$$\widehat{S}_{F,ad}(E_n,E)=\widehat{S}_{F,ad}(E,E_n)=\widehat{S}_{0,n}(E).$$
(16)
The adiabatic Floquet scattering matrix consists thus for each incident energy of a block of dimension of $`(2N_rn_{max})^2`$ where $`N_r`$ is the number of leads and $`n_{max}`$ is the maximum number of side bands needed for an accurate description of the quantities of interest.

The above equation allows us to formulate the following adiabaticity criterion:. If the Floquet scattering matrix $`\widehat{S}_F`$ changes only a little when the energy $`E`$ changes by $`n_{max}\mathrm{}\omega `$ then the adiabatic approximation can be applied and $`\widehat{S}_F\widehat{S}_{F,ad}`$. Note that because of Eq.(16) the same criterion can be applied to the scattering matrix $`\widehat{S}_0(E)`$.

Let us now calculate the directed charge current $`I_{ad}`$ and a heat flow $`I_{E,ad}`$ in the adiabatic limit. Substituting Eq.(16) into Eq.(8) we find the current flowing in lead $`\alpha `$ under an adiabatic change of parameters (we put $`f^{(in)}=f_0`$):

$$\begin{array}{c}I_{ad,\alpha }=\frac{e\omega }{2\pi }\underset{0}{\overset{\mathrm{}}{}}𝑑E\left(\frac{f_0(E)}{E}\right)\\ \\ \times \underset{\beta }{}\underset{n}{}n|S_{0,\alpha \beta ,n}(E)|^2.\end{array}$$
(17)
We see that the current $`I_{ad}`$ is due to photon-assisted processes (the current depends on the intensity of side bands which is proportional to $`|S_{0,\alpha \beta ,n}(E)|^2`$). Therefore even in the ”adiabatic” limit pumping is, strictly speaking, a nonadiabatic phenomenon.

Note that formally Eq. (17) is obtained for $`\mathrm{}\omega k_BT`$. However it gives the correct answer (in the adiabatic limit) for an arbitrary ratio of the frequency $`\omega `$ to the temperature $`T`$.

Using Eq.(14) we can rewrite Eq.(17) as follows

$$\begin{array}{c}I_{ad,\alpha }=i\frac{e\omega }{4\pi ^2}\underset{0}{\overset{2\pi /\omega }{}}𝑑t\underset{0}{\overset{\mathrm{}}{}}𝑑E\left(\frac{f_0(E)}{E}\right)\\ \\ \times \left(\frac{\widehat{S}_0(E,t)}{t}\widehat{S}_0^{}(E,t)\right)_{\alpha \alpha }.\end{array}$$
(18)
This is Brouwer’s formula . We emphasize that Eq.(17) and Eq.(18) express a (nonadiabatic in nature) pumped current in terms of the scattering matrix $`\widehat{S}_0`$ for a ”frozen” scatterer with time-dependent parameters (see Eq.(13)). This is correct if the scattering matrix $`\widehat{S}_0`$ can be taken to be energy independent on the scale of the order of $`n_{max}\mathrm{}\omega `$. The last conclusion is correct irrespective of the amplitude of the oscillating parameters. But $`n_{max}`$ might very well depend on whether one is in the weak amplitude limit or in the strong amplitude limit.

Next consider the heat flow $`I_{E,ad,\alpha }`$ produced by the adiabatically oscillating scatterer in the lead $`\alpha `$. Substituting the adiabatic scattering matrix Eq.(16). into Eq.(12) we find

$$\begin{array}{c}I_{E,ad,\alpha }=\frac{1}{h}\underset{\mathrm{}}{\overset{\mathrm{}}{}}𝑑ϵ\underset{\beta }{}\underset{n}{}|S_{0,\alpha \beta ,n}(\mu +ϵ)|^2\\ \\ \times \frac{\left(ϵ+\frac{n\mathrm{}\omega }{2}\right)\mathrm{sinh}\left(\frac{n\mathrm{}\omega }{2k_BT}\right)}{\mathrm{cosh}\left(\frac{ϵ}{k_BT}\right)+\mathrm{cosh}\left(\frac{n\mathrm{}\omega }{2k_BT}\right)}.\end{array}$$
(19)
At zero temperature (or if the temperature $`T`$ is less than the relevant energy scale for the scattering matrix $`\widehat{S}_0`$) the above equation can be greatly simplified. In this case the heat flow can be expressed in terms of an on-shell, time-dependent adiabatic scattering matrix $`\widehat{S}_0(\mu ,t)`$:

$$\begin{array}{c}I_{E,ad,\alpha }=\frac{\mathrm{}\omega }{8\pi ^2}\underset{0}{\overset{2\pi /\omega }{}}𝑑t\hfill \\ \\ \times \left(\frac{^2\widehat{S}_0(\mu ,t)}{t^2}\widehat{S}_0^{}(\mu ,t)\right)_{\alpha \alpha }.\hfill \end{array}$$
(20)
This expression coincides with that given by Avron et al. .

In the adiabatic limit it is possible to express the quantities of interest in terms of the stationary scattering matrix. This permits the very general and useful expressions discussed in this paragraph. In contrast, in the non-adiabatic case, the information contained in the stationary scattering matrix is by definition not sufficient. As a consequence, in order to address the non-adiabatic case, we have now to consider a very specific example.

## IV Charge and heat flows produced by an oscillating double barrier

To apply the general formalism of Sec.II beyond the adiabatic approximation we investigate a simple model . It is a one-dimensional scatterer consisting of two delta function barriers oscillating with frequency $`\omega `$ and located at $`x=L/2`$ and $`x=L/2`$.

To calculate the Floquet matrix we consider scattering of electrons with energy $`E`$ coming from $`x=\mathrm{}`$. The time-dependent Schrödinger equation for an electron wave function $`\mathrm{\Psi }(x,t)`$ reads

$$\begin{array}{c}i\mathrm{}\frac{\mathrm{\Psi }(x,t)}{t}=\widehat{H}(x,t)\mathrm{\Psi }(x,t),\hfill \\ \\ \widehat{H}(x,t)=\frac{\mathrm{}^2}{2m}\frac{^2}{x^2}+V(x,t),\hfill \\ \\ V(x,t)=V_1(t)\delta (x+\frac{L}{2})+V_2(t)\delta (x\frac{L}{2}),\hfill \\ \\ V_i(t)=V_{0i}+2V_{1i}\mathrm{cos}(\omega t+\phi _i),i=1,2.\hfill \end{array}$$
(21)
Appendix A gives the exact solution of this model. We now use this solution to calculate the charge current flowing through (as well as the heat current flowing from) the oscillating double barrier potential (see Eq.(21)) connecting two reservoirs with the same temperature $`T=0`$ and electrochemical potential $`\mu `$. We use the units $`2m=\mathrm{}=e=1`$. To be definite we consider the current $`I_1`$ ($`I_{E,1}`$) flowing to the left ($`x\mathrm{}`$). The currents are calculated, using Eq.(8) and Eq.(12) with the Floquet scattering matrix determined as shown in Appendix A. We compare these currents with the adiabatic currents $`I_{ad,1}`$ and $`I_{E,ad,1}`$ using the Brouwer formula Eq.(18) and its analog Eq.(20).

For reference we write down the scattering matrix $`\widehat{S}_0(E)`$ which we need to calculate the adiabatic currents:

$$\widehat{S}_0=\frac{e^{ikL}}{\mathrm{\Delta }}\left(\begin{array}{cc}\xi +2\frac{p_2}{k}\mathrm{sin}(kL)& 1\\ & \\ 1& \xi +2\frac{p_1}{k}\mathrm{sin}(kL)\end{array}\right).$$
(22)
Here $`k=\sqrt{\frac{2m}{\mathrm{}^2}E}`$; $`p_j=V_jm/\mathrm{}^2`$ (j = 1,2); $`\xi `$ = $`(1\mathrm{\Delta })e^{ikL}`$; $`\mathrm{\Delta }`$ = $`1+\frac{p_1p_2}{k^2}(e^{2ikL}1)`$ \+ $`i\frac{p_1+p_2}{k}`$. Of interest is a comparison of the exact Floquet scattering matrix and adiabatic theory in the limit of small and large frequencies. Furthermore, it is interesting to compare the symmetry conditions for pumping for the two theories.

### A Adiabatic limit: $`\omega 0`$

¿From our calculations it follows that in the limit of $`\omega 0`$ the adiabatic current $`I_{ad}`$ gives a good approximation for the pumped current irrespective of the amplitude of the oscillating potentials $`V_{1i}`$. However the criterion of adiabaticity depends strongly on the ratio $`V_{1i}/\mathrm{}\omega `$.

As we have seen, the adiabatic approximation is valid if the scattering matrix $`\widehat{S}_0`$ is energy independent on the scale of the order of $`n_{max}\mathrm{}\omega `$. Here $`n_{max}`$ (specified below) is the number of side bands with noticeable amplitude excited by the oscillating scatterer. Let us denote by $`\delta E`$ the energy scale over which the scattering matrix $`\widehat{S}_0`$ changes significantly. For example, close to a resonance $`\delta E`$ is of the order of the width $`\delta `$ of a resonance. Then the adiabatic approximation is valid if $`n_{max}\mathrm{}\omega \delta E`$. Hence the larger the number $`n_{max}`$ of excited side bands the smaller the frequencies for which the adiabatic approximation is valid.

As it is well known  the number $`n_{max}`$ of excited side bands depends on the ratio of the amplitude of the oscillating potential $`V_{1i}`$ to the frequency $`\omega `$. In the small amplitude limit $`V_{1i}0`$ only the first side bands are excited: $`n_{max}=1`$. But for strong pumping the number $`n_{max}`$ is large and the adiabatic approximation is valid only at smaller frequencies. As a consequence for a given finite frequency $`\omega `$ the deviation of the actual pumped current $`I(\omega )`$ from the adiabatic one $`I_{ad}`$ increases with increasing pumping amplitude. This fact is illustrated in Fig.1.

Before concluding this subsection we would like to emphasize the following. For strong pumping the adiabatic approximation is still valid at sufficiently small frequencies despite the excitation of a large number of side bands. The Fourier harmonics of the scattering matrix $`\widehat{S}_0`$ define the amplitudes of side bands for strong (Eq.(16)) as well as for weak  adiabatic pumping. Thus the scattering matrix $`\widehat{S}_0`$ completely defines the kinetics and, in particular, the quantum statistical correlation properties of an adiabatic pump. The current noise of a weak (small amplitude) pump was considered in Ref.. The noise in a strong amplitude adiabatic pump will be presented elsewhere.

### B Large frequency limit

At large frequencies $`\omega >\delta `$ (where $`\delta `$ is the width of a resonance) the pumped current $`I`$ differs considerably from the adiabatic one $`I_{ad}`$.

It is known that for a double barrier structure at $`\omega 0`$ the pumped dc current $`I_{ad}`$  and the generated heat flow $`I_{E,ad}`$  show a resonance-like behavior as a function of the Fermi energy $`\mu `$ (see Fig.2). The pumped current and heat peak when the Fermi energy is close to the transmission resonance. This is because in the adiabatic limit only the particles close to the Fermi level contribute to the charge (Eq.(18)) and energy (Eq.(20)) transfer.

With increasing pumped frequency $`\omega \delta `$ the particles within a wider energy interval come into play. The dependence on the Fermi energy is smoothed away. However a resonance-like dependence on the frequency $`\omega `$ arises.

In Fig.3 we depict the dependence of the charge $`\delta Q`$ = $`(2\pi )I/\omega `$ pumped during a cycle on the frequency $`\omega `$. The pumped charge peaks when the energy quantum $`\mathrm{}\omega `$ equals one (or several) level spacings $`\mathrm{\Delta }`$ of a double barrier structure. In the example used for numerical calculations the level spacing $`\mathrm{\Delta }`$ near the Fermi energy is constant with good accuracy. The smaller peaks correspond to many ($`n=2,3,\mathrm{}`$) photon processes for which $`n\mathrm{}\omega =m\mathrm{\Delta }`$ (m = 1,2,…). We see that the single photon processes dominate over the many photon processes.

At arbitrary ($`n\mathrm{}\omega m\mathrm{\Delta }`$) but large ($`\mathrm{}\omega \delta `$) frequencies, due to interference (inside the double barrier structure) only the main component ($`n=0`$) of the Floquet state Eq.(23) has a significant amplitude between the barriers. This component corresponds to the eigenfunction of the time independent problem (with barriers $`V_{01}`$ and $`V_{02}`$). The side bands ($`n0`$) do not participate in the transmission through the system (more precisely, their contribution is small). As a result the pumped current is greatly reduced compared to the adiabatic case. It should be noted that in the adiabatic case the side bands do contribute to the transmission since for $`\omega \delta `$ they lie at the same transmission resonance as the main component.

On the other hand at some particular values of a frequency $`\mathrm{}\omega =m\mathrm{\Delta }`$ the substates $`\psi _{\pm k}`$ ($`k=1,2,3,\mathrm{}`$) of the Floquet state Eq.(23) become large and additional channels for the transmission through the system open up. This leads to the increase of the pumped charge (see main peaks in Fig.3).

Interestingly, as shown in Fig.3), the pumped current reverses sign as a function of frequency. These sign reversals can be understood in the following way: If two potential barriers oscillate then the pump effect arises as a consequence of an interference between two amplitudes . The first amplitude $`𝒜_1`$ corresponds to particles which propagate through the double barrier structure and absorb (emit) the energy quantum $`\mathrm{}\omega `$ (or several quanta $`n\mathrm{}\omega `$) at the first barrier. The second amplitude $`𝒜_2`$ corresponds to the same propagation with absorption (emission) at the second barrier. As it is known  the adiabatically pumped current is an odd function of the phase difference of the oscillating potentials $`\mathrm{\Delta }\phi `$ = $`\phi _1`$ \- $`\phi _2`$ (see Eq.(21)). In fact this is the phase difference of the two amplitudes $`𝒜_1`$ and $`𝒜_2`$. In the nonadiabatic case there exists an additional contribution coming from the spatial phase difference

$$\mathrm{\Delta }\phi _x(k_1k_2)L.$$

Thus at large frequencies the pumped current depends on the sum $`\mathrm{\Delta }\phi `$ \+ $`\mathrm{\Delta }\phi _x`$.

To clarify the appearance of $`\mathrm{\Delta }\phi _x`$ let us consider particles with energy $`E`$ going from the left to the right and absorbing the quantum $`\mathrm{}\omega `$. The amplitude $`𝒜_1`$ corresponds to particles which absorb an energy $`\mathrm{}\omega `$ close to the left barrier ($`V_1`$) and traverse the system at energy $`E+\mathrm{}\omega `$. Hence $`k_1\sqrt{E+\mathrm{}\omega }`$. On the other hand the amplitude $`𝒜_2`$ corresponds to particles which absorb an energy $`\mathrm{}\omega `$ close to the right barrier ($`V_2`$) and thus traverse the system at energy $`E`$. They thus propagate with wave vector $`k_2\sqrt{E}`$. Because $`k_1k_2`$ there is a phase difference $`\mathrm{\Delta }\phi _x`$ between the amplitudes $`𝒜_1`$ and $`𝒜_2`$. Note that in the adiabatic case (i.e., at $`\omega 0`$) $`\mathrm{\Delta }\phi _x`$ vanishes.

Now we show that in the case under consideration the phase difference $`\mathrm{\Delta }\phi _x`$ is at the origin of the sign reversal of the pumped current at consecutive peaks (see Fig.3). Close to the $`m^{th}`$ main peak of $`I(\omega )`$ the frequency is determined by $`\mathrm{}\omega m\mathrm{\Delta }`$. Taking into account that only the particles with energy close to the resonance contribute to the transmission and mainly single photon processes are important we can estimate $`\mathrm{\Delta }\phi _x`$ as follows. If the amplitude $`𝒜_1`$ corresponds to the propagation through some resonance (say $`E_{l+m}`$) then $`𝒜_2`$ corresponds to the propagation through the resonance $`E_l`$. At the $`l^{th}`$ transmission resonance it is $`k^{(l)}L=l\pi `$ the additional phase difference is $`\mathrm{\Delta }\phi _x=m\pi `$. Thus the sum $`\mathrm{\Delta }\phi `$ \+ $`\mathrm{\Delta }\phi _x`$ changes by $`\pi `$ as we pass from one peak to another ($`\mathrm{\Delta }m=1`$) thus giving rise to the sign reversal .

In addition interference inside the double barrier structure manifests itself as a resonance-like dependence of the heat production (see Fig.4). We can see that $`I_E(\omega )`$ contains additional peaks at $`\mathrm{}\omega `$ = $`(m+\frac{1}{2})\mathrm{\Delta }`$. This fact emphasizes the striking difference between the processes which are responsible for charge pumping and for heat production .

### C Quantum pumping and broken symmetry

Our results show that the symmetry of the scattering problem is important for the pump effect.

The quantum pump effect arises due to the different, interfering excitation histories of carriers traversing the sample . To ”extract” directed currents from the uniform environment the left-right symmetry (LRS) for carriers traversing the sample needs to be broken. There are, at least, two ways to break the LRS. The first one is breaking the spatial symmetry (SS). The second one is breaking the time reversal symmetry (TRS).

We restrict our considerations to systems which in the absence of the time-dependent perturbations needed for pumping are in an equilibrium state. In particular this means that without the presence of magnetic fields (etc.) the system is time reversal invariant. Note that the time independent scattering problem is insensitive to the presence (absence) of the spatial symmetry of the scatterer (in the sense that the transmission probability is invariant under the spatial inversion $`xx`$ irrespectively of the SS of the scattering potential).

Interestingly, the symmetry conditions for pumping depend on the frequency. In the adiabatic limit $`\omega 0`$ (see Sec.III) the scattering problem is fully characterized by the scattering matrix $`\widehat{S}_0`$. In contrast to the conductance, adiabatic pumping is sensitive to the symmetry of the scattering matrix. To have an adiabatic pumping effect we need to break both the spatial symmetry and the time-reversal symmetry.

It is important to distinguish the symmetry of the equilibrium problem and the symmetry of the full problem. In the adiabatic case, the symmetry of the time-independent problem is irrelevant: What matters, is that the system in presence of perturbations breaks both the spatial symmetry and the time reversal symmetry. In many examples (for instance the two barrier problem considered above) it is not possible to break the time-reversal symmetry without at the same time breaking the spatial symmetry. However, examples which are spatially symmetric and have broken time-reversal invariance can be constructed. For instance, we could consider the problem of two barriers oscillating in synchronism and as a second perturbation which oscillates with a phase lag choose the potential between the two barriers. This would be an example of a two parameter problem which does not generate a pumped current. The important point we would like to emphasize is that the symmetry conditions for the non-adiabatic pump are different.

In the strongly nonadiabatic limit the scattering problem is described by the Floquet scattering matrix $`\widehat{S}_F(E^{(out)},E^{(in)})`$ dependent on both the incident $`E^{(in)}`$ and outgoing $`E^{(out)}`$ energies. In general, the Floquet scattering matrix is sensitive to the spatial symmetry of the scatterer. If the scattering problem is without SS (i.e., when the scattering potential is not invariant under the change $`xx`$) then

$$\begin{array}{c}|S_{F,\alpha \beta }(E^{(out)},E^{(in)})|^2\\ \\ |S_{F,\beta \alpha }(E^{(out)},E^{(in)})|^2,\end{array}$$

and the oscillating scatterer can pump a dc current (see Eq.(9)). Thus to obtain a pump effect in the strongly nonadiabatic limit there are, at least, two possibilities. Either the spatial symmetry or the time reversal symmetry (or both SS and TRS together) have to be broken. The role of spatial symmetry breaking for the nonadiabatic pump effect was emphasized in Ref..

The above reasoning is illustrated in Fig.5 where the dependence of a pumped current on the frequency is shown for three generic situations. First, two potentials oscillate with the same amplitude but out of phase (solid line). In this case the perturbation breaks both the spatial symmetry and the time reversal symmetry. Second, only one potential oscillates (dotted line). The TRS is present but the SS is dynamically broken. The system is able to pump a dc current like in the previous case. Third, two potentials oscillate in phase and with the same amplitudes (dash-dotted line). Now both the TRS and SS are preserved and the system under consideration does not exhibit a directed current.

## V Discussion and conclusion

In this work we have developed the Floquet scattering matrix approach to parametric pumping in phase coherent mesoscopic systems of noninteracting electrons. We have calculated the distribution function, the dc current and the heat flow at arbitrary pumping amplitudes and frequency.

The Floquet scattering matrix describes naturally the existence of the side bands $`E_n=E+n\mathrm{}\omega `$ of particles leaving the pump. These side bands correspond to nonequilibrium particles generated by the pump which carry the heat from the oscillating scatterer to the reservoirs and (under some conditions) transfer charge between the reservoirs . If the pumping amplitude is small only the first ($`n=\pm 1`$) side bands are excited. But for strong pumping the number of excited side bands is large: $`n_{max}1`$. This number together with the frequency defines the energy scale $`n_{max}\mathrm{}\omega `$ characteristic of the pump problem. In particular this is important for the analysis of the conditions under which the adiabatic approach to pumping is valid. If the scattering matrix is energy independent on the scale of the order of $`n_{max}\mathrm{}\omega `$ then the adiabatic approximation can be applied. In this case the elements of the Floquet scattering matrix are given by the corresponding Fourier components of the stationary scattering matrix $`\widehat{S}_0`$ (see Eq.(16)) with parameters taken to be dependent on time.

The existence of the pump effect is directly related to the symmetry of the scattering problem. In the adiabatic case ($`\omega 0`$) only a scatterer without spatial and time reversal symmetry can produce a directed current. This conclusion applies irrespectively of the (stationary) spatial symmetry of the scatterer. On the other hand in the nonadiabatic case (at large pumping frequency) to achieve pumping we need a scatterer with either broken spatial or time reversal symmetry. Hence at large operating frequencies even the scatterer with a single oscillating parameter can show a pump effect if only the spatial symmetry is broken .

To emphasize the main physics underlying pumping we have considered an exactly solvable model. Namely a one-dimensional scatterer consisting of two oscillating delta-function barriers Eq.(21) separated by a distance $`L`$.

If two barriers oscillate out of phase (broken TRS) the basic process leading to the pump effect is an interference of two quantum mechanical amplitudes. These amplitudes $`𝒜_1`$ and $`𝒜_2`$ corresponds to particles which traverse the scatterer in the same direction but gain (or lose) a modulation energy in the vicinity either of the first or of the second barrier, respectively . This becomes clearer at large frequencies $`\omega `$ because of the following. When the energy quantum $`\mathrm{}\omega `$ equals one (or several, say $`m`$) level spacing $`\mathrm{\Delta }`$ of a double barrier structure the pumped current peaks (see Fig.3). The particles with energy close to the resonance of a double barrier structure give the main contribution to the pump current. Thus the amplitudes $`𝒜_1`$ and $`𝒜_2`$ correspond to traversal of the system through different (say, $`(l+m)^{th}`$ and $`l^{th}`$) resonances. This gives rise to an additional phase difference $`\mathrm{\Delta }\phi _x`$ = $`(k_1k_2)L`$ between the amplitudes $`𝒜_1`$ and $`𝒜_2`$. For the resonant double barrier considered here the phase difference between consecutive resonances is $`\pi `$ and thus $`\mathrm{\Delta }\phi _x=m\pi `$. As a result the consecutive peaks of the pumped current have opposite sign (see Fig.3).

We conclude that if the time reversal symmetry is broken (i.e., $`\mathrm{\Delta }\phi 0`$), the nonadiabatic pumped current contains direct information on the phase of the transmission coefficient through the mesoscopic sample. We emphasize that here we have a single-connected geometry in contrast with the case where the mesoscopic sample (a quantum dot) is embedded in one arm of the ring .

We have also considered the effect of a spatial asymmetry. If only one potential (say, $`V_1`$) oscillates then the spatial symmetry is (dynamically) broken. In this case the (nonadiabatic) pump effect arises because the probability for exciting side bands $`|S_{F,\alpha \beta }(E_n,E)|^2`$ depends on the direction (i.e., on the order of indexes $`\alpha `$ and $`\beta `$) (see Eq.(9)). The particles going from the left to the right first absorb (emit) the energy $`n\mathrm{}\omega `$ = $`E_nE`$ and then pass through the double barrier structure (at energy $`E_n`$). Thus the corresponding probability is

$$|S_{F,21}(E_n,E)|^2T(E_n).$$

Here $`T(E)`$ is a transmission probability. On the other hand the particles going from the right to the left pass the system before they absorb (or emit) the energy in the vicinity of the oscillating barrier $`1`$. Hence we have

$$|S_{F,12}(E_n,E)|^2T(E).$$

Thus the pumped current Eq.(9) depends on the differences $`T(E)T(E\pm n\mathrm{}\omega )`$.

Note that for a system with equidistant spectrum (that is the case under consideration) the transmission probability $`T(E)`$ is periodic in energy $`E`$ with the period of $`\mathrm{\Delta }`$. Hence if $`\mathrm{}\omega =m\mathrm{\Delta }`$ then it is $`T(E)T(E\pm \mathrm{}\omega )`$ and the current is close to zero in accordance with Fig.5 (a dotted curve; $`m=3`$).

We have presented a scattering theory of quantum pumping based on the Floquet theorem. This approach allows a description of the kinetics of both adiabatic and nonadiabatic quantum pumps. In particular the quantum statistical correlation properties (noise) can be considered in analogy with the weak amplitude adiabatic quantum pump . The investigation of noise is especially important in the context of quantized (i.e., noiseless) charge pumping  and will be presented elsewhere.

ACKNOWLEDGEMENT

This work is supported by the Swiss National Science Foundation.

## VI Appendix A: Oscillating double barrier. An exact solution

In this Appendix we determine the Floquet scattering matrix of the Schrödinger equation (21) of two oscillating delta function potentials located at $`x=L/2`$ and $`x=L/2`$. To find the Floquet matrix we consider scattering of electrons with energy $`E`$ coming from $`x=\mathrm{}`$. To solve (21) we use the Floquet functions method . Since the Hamiltonian $`\widehat{H}`$ depends on time, the system has no stationary eigenstates. However, the Floquet theorem tells us that because the Hamiltonian is periodic in time the eigenstates of Eq.(21) can be represented as a superposition of wave functions with energies shifted by $`n\mathrm{}\omega `$:

$$\mathrm{\Psi }_E(x,t)=e^{iEt/\mathrm{}}\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}\psi _n(x)e^{in\omega t}.$$
(23)
Away from the points $`x=L/2`$ and $`x=L/2`$ the functions $`\psi _n(x)`$ are a superposition of plane waves

$$\psi _n(x)=a_ne^{ik_nx}+b_ne^{ik_nx},$$
(24)
where

$$k_n=\sqrt{\frac{2m}{\mathrm{}^2}\left(E+n\mathrm{}\omega \right)}$$

with $`Re[k_n]0`$ and $`Im[k_n]0`$. At $`x=L/2`$ and $`x=L/2`$ we have the boundary conditions

$$\begin{array}{c}\mathrm{\Psi }_E(\frac{L}{2}0,t)=\mathrm{\Psi }_E(\frac{L}{2}+0,t),\\ \\ \mathrm{\Psi }_E(\frac{L}{2}0,t)=\mathrm{\Psi }_E(\frac{L}{2}+0,t),\\ \\ \frac{\mathrm{\Psi }_E(x,t)}{x}|_{x=\frac{L}{2}+0}\frac{\mathrm{\Psi }_E(x,t)}{x}|_{x=\frac{L}{2}0}\\ \\ =\frac{2m}{\mathrm{}^2}V_{01}(t)\mathrm{\Psi }_E(\frac{L}{2},t),\\ \\ \frac{\mathrm{\Psi }_E(x,t)}{x}|_{x=\frac{L}{2}+0}\frac{\mathrm{\Psi }_E(x,t)}{x}|_{x=\frac{L}{2}0}\\ \\ =\frac{2m}{\mathrm{}^2}V_{02}(t)\mathrm{\Psi }_E(\frac{L}{2},t).\end{array}$$
(25)
In the case under consideration (i.e., when an electron with energy $`E`$ comes from $`x=\mathrm{}`$) the functions $`\psi _n(x)`$ are

$$\begin{array}{c}\psi _n(x<0)=\delta _{n,0}e^{ik_nx}+R_ne^{ik_nx},\hfill \\ \\ \psi _n(0<x<L)=A_ne^{ik_nx}+B_ne^{ik_nx},\hfill \\ \\ \psi _n(x>L)=T_ne^{ik_nx}.\hfill \end{array}$$
(26)
Here the coefficients $`R_n`$ and $`T_n`$ for propagating modes (for which $`E_n>0`$) are the amplitudes of reflection from or transmission through the double barrier system absorbing ($`n>0`$) or emitting ($`n<0`$) an energy $`|n|\mathrm{}\omega `$, respectively.

Substituting Eq.(23) and Eq.(26) into the boundary conditions Eq.(25) we obtain

$$\begin{array}{c}R_n=\left(A_n\delta _{n,0}\right)e^{ik_nL}+B_n,\hfill \\ \\ T_n=A_n+B_ne^{ik_nL},\hfill \end{array}$$
(27)
where the coefficients $`A_n`$ and $`B_n`$ are subject to recursive equations. It is convenient to represent these equations in a matrix form

$$\begin{array}{c}\widehat{U}_n\left(\begin{array}{c}A_n\hfill \\ B_n\hfill \end{array}\right)+\delta _{n,0}\left(\begin{array}{c}ik_0e^{ik_0\frac{L}{2}}\\ 0\end{array}\right)\\ \\ =\widehat{D}_{n+1}\left(\begin{array}{c}A_{n+1}\hfill \\ B_{n+1}\hfill \end{array}\right)+\widehat{I}_{n1}\left(\begin{array}{c}A_{n1}\hfill \\ B_{n1}\hfill \end{array}\right).\end{array}$$
(28)
Here we have introduced the matrices

$$\widehat{U}_n=\left(\begin{array}{cc}(ik_np_{01})e^{ik_n\frac{L}{2}}\hfill & \hfill p_{01}e^{ik_n\frac{L}{2}}\\ & \\ p_{02}e^{ik_n\frac{L}{2}}\hfill & \hfill (ik_np_{02})e^{ik_n\frac{L}{2}}\end{array}\right),$$
(29)
$$\widehat{D}_{n+1}=\left(\begin{array}{cc}p_{11}e^{i\phi _1}e^{ik_{n+1}\frac{L}{2}}\hfill & \hfill p_{11}e^{i\phi _1}e^{ik_{n+1}\frac{L}{2}}\\ & \\ p_{12}e^{i\phi _2}e^{ik_{n+1}\frac{L}{2}}\hfill & \hfill p_{12}e^{i\phi _2}e^{ik_{n+1}\frac{L}{2}}\end{array}\right),$$
(30)
$$\widehat{I}_{n1}=\left(\begin{array}{cc}p_{11}e^{i\phi _1}e^{ik_{n1}\frac{L}{2}}\hfill & \hfill p_{11}e^{i\phi _1}e^{ik_{n1}\frac{L}{2}}\\ & \\ p_{12}e^{i\phi _2}e^{ik_{n1}\frac{L}{2}}\hfill & \hfill p_{12}e^{i\phi _2}e^{ik_{n1}\frac{L}{2}}\end{array}\right),$$
(31)
where the parameters are: $`p_{ij}=V_{ij}m/\mathrm{}^2`$, $`i=0,1`$, $`j=1,2`$.

To solve the equation (28) we have generalized the method used in Ref. for a single oscillating delta-function potential.

First of all we consider Eq.(28) at $`n>0`$. We suppose that there exist matrices $`\widehat{X}_n`$ such that

$$\left(\begin{array}{c}A_n\hfill \\ B_n\hfill \end{array}\right)=\widehat{X}_n\left(\begin{array}{c}A_{n1}\hfill \\ B_{n1}\hfill \end{array}\right).$$
(32)
Then substituting the above constraint into Eq.(28) we obtain a simple recursive equation for the matrices $`\widehat{X}_n`$ ($`n>0`$)

$$\widehat{X}_n=\left(\widehat{U}_n\widehat{D}_{n+1}\widehat{X}_{n+1}\right)^1\widehat{I}_{n1}.$$
(33)
Using these matrices we can express all the coefficients $`A_n`$, $`B_n`$ ($`n>0`$) in terms of $`A_0`$ and $`B_0`$ only

$$\left(\begin{array}{c}A_n\hfill \\ B_n\hfill \end{array}\right)=\widehat{X}_n\mathrm{}\widehat{X}_1\left(\begin{array}{c}A_0\hfill \\ B_0\hfill \end{array}\right),n>0.$$
(34)
Further we consider $`n<0`$ and introduce the matrices $`\widehat{Y}_n`$

$$\left(\begin{array}{c}A_n\hfill \\ B_n\hfill \end{array}\right)=\widehat{Y}_n\left(\begin{array}{c}A_{n+1}\hfill \\ B_{n+1}\hfill \end{array}\right).$$
(35)
The corresponding recursive equation for the matrices $`\widehat{Y}_n`$ is

$$\widehat{Y}_n=\left(\widehat{U}_n\widehat{I}_{n1}\widehat{Y}_{n1}\right)^1\widehat{D}_{n+1}.$$
(36)
Using these matrices we express the coefficients $`A_n`$, $`B_n`$ ($`n<0`$) in terms of $`A_0`$ and $`B_0`$:

$$\left(\begin{array}{c}A_n\hfill \\ B_n\hfill \end{array}\right)=\widehat{Y}_n\mathrm{}\widehat{Y}_1\left(\begin{array}{c}A_0\hfill \\ B_0\hfill \end{array}\right),n<0.$$
(37)
As a final step we consider Eq.(28) at $`n=0`$. After simple manipulations we find the coefficients $`A_0`$ and $`B_0`$:

$$\begin{array}{c}\left(\begin{array}{c}A_0\hfill \\ B_0\hfill \end{array}\right)=\left(\widehat{U}_0\widehat{D}_1\widehat{X}_1\widehat{I}_1\widehat{Y}_1\right)^1\\ \\ \times \left(\begin{array}{c}ik_0e^{ik_0\frac{L}{2}}\\ 0\end{array}\right).\end{array}$$
(38)
Thus, using the solutions of the recurrent equations Eq.(33) and Eq.(36) we can calculate all the coefficients $`A_n`$, $`B_n`$ (see Eq.(34), Eq.(37), and Eq.(38)). Note that in each particular case we need to take into account only the limited number $`|n|<n_{max}`$ of side bands and thus we can put $`\widehat{X}_{n_{max}+1}0`$ and $`\widehat{Y}_{n_{max}1}0`$.

The coefficients $`R_n`$ and $`T_n`$ can be calculated with the help of Eq.(27). For the propagating modes ($`E_nE+n\mathrm{}\omega >0`$) of interest here, these coefficients determine the elements of the Floquet scattering matrix

$$\begin{array}{c}|S_{F,11}(E_n,E)|^2=\frac{k_n}{k_0}|R_n|^2,\\ \\ |S_{F,21}(E_n,E)|^2=\frac{k_n}{k_0}|T_n|^2.\end{array}$$
(39)
Here the indexes $`1`$ and $`2`$ correspond to the left and right reservoirs, respectively. To obtain the matrix elements $`S_{F,22}`$ and $`S_{F,12}`$ we need to solve the same problem with plane waves coming from the right.


format: text

text:
# Axisymmetric versus Non-axisymmetric Vortices in Spinor Bose-Einstein Condensates

## I Introduction

The experimental achievement of Bose-Einstein condensation (BEC) in the trapped atomic cloudscornell ; bradley ; davis  has opened up a novel field to investigate fundamental problems such as the relationship between the superfluidity and BEC. Owing to recent advances of experimental techniques, several groups have succeeded in creating quantized vortices with various procedures in the magnetically trapped BECmatthews ; madison ; abo ; haljan ; hodby ; leanhardt , where the condensate is described by a scalar order parameter. Furthermore, the atomic gases with the hyperfine spin $`F=1`$ called “spinor BEC” have been Bose-condensed via optical methodsstenger ; barrett  which can keep the atomic “spin” states degenerate and active. As shown recently by Klausen et al.klausen , the spin-dependent interaction of two 87Rb atoms is ferromagnetic. Thus, we now have spinor BEC’s with both antiferromagnetic (23Na)stenger  and ferromagnetic interactions.

Such scalar and spinor BEC’s are analogous to superfluid 3He and 4He. These superfluid Heliums, however, have rather strong interactions. Indeed, the condensate fraction in superfluid 4He is only 10$`\%`$ of the total. By contrast, BEC of the atomic gases have advantages for both theoretical and experimental treatments due to their weak interactions; here almost all the atoms are able to be Bose-condensed. It is possible to directly observe dynamical behaviors of the condensate with optical methods, providing us an opportunity to quantitatively investigate the new quantum fluiddalfovo .

The standard Hamiltonian for the spinor BEC have been introduced by Ohmi and Machidaohmi , and Hoho , who pointed out the richness of the exotic topological defects. Topological structures, such as skyrmion, meron, Mermin-Ho (Anderson-Toulouse) texture, and monopole, play an important role in various fields of physics. They provide a common framework to connect diverse field, thereby enhancing mutual understandingjayaraman ; leggett . Specially, since $`F=1`$ ferromagnetic BEC can be described by order parameters similar to the superfluid 3He-A, the coreless Mermin-Ho vortex may be favored in the ferromagnetic BECmermin ; anderson .

Similar topological structures, called skyrmion in general, have been proposed in the spinor BEC. Al Khawaja and Stoof stoof  studied a skyrmion in the $`F=1`$ ferromagnetic BEC and concluded that it is not a thermodynamically stable object without rotation. By considering the effect of the external rotation, however, we have shown recentlymizushima  that this topological defect can be stable. Ingenious proposals have been mademarzlin ; anglin ; kurihara ; martikainen  on how to create it and detect it. Yipyip  has performed a systematic study on vortex structures and presented several axisymmetric and non-axisymmetric vortices for $`F=1`$ antiferromagnetic BEC. Recently, Isoshima et al.isoshima1 ; isoshima2  have carried out an extensive study of axisymmetric vortices to provide a vortex phase diagram in a plane of the rotation and the magnetization for both the antiferromagnetic and ferromagnetic cases.

In this paper we examine the stability of various vortices for both the $`F=1`$ antiferromagnetic and ferromagnetic BEC trapped in a two-dimensional harmonic potential. We have removed the previously imposed restriction in the axisymmetric case that winding numbers are less than or equal to unity. The continuous vortices such as the Mermin-Ho vortex will also be shown to be favored over the singular onesisoshima1 ; isoshima2  and the other non-axisymmetric onesyip . We demonstrate the stability of such vortices and discuss differences between the axisymmetric and non-axisymmetric configurations. We also determine the vortex phase diagrams for the antiferromagnetic and ferromagnetic cases. By comparing the relative free energies of the possible vortex configurations and the phase-separated state which may occur in the ferromagnetic situation, the validity of assuming uniformity along $`z`$ direction will be checked.

This paper is organized as follows. In Sec. II, we first present the extended Gross-Pitaevskii equation for the spinor BEC, and then explain the numerical procedure to find local minima of the energy functional. Section III enumerates possible vortices for axisymmetric and non-axisymmetric cases, and then discuss the stability of the Mermin-Ho vortex. Section IV presents the phase diagram for the ferromagnetic and antiferromagnetic cases in the plane of external rotation $`\mathrm{\Omega }`$ vs total magnetization $`M`$ obtained by comparing the free energies. Here, the qualitative differences between the axisymmetric and non-axisymmetric vortices are discussed by showing the $`\mathrm{\Omega }`$\- and $`M`$-dependence. The final section is devoted to conclusions and discussions.

## II Theoretical Formulation

### II.1 Extended Gross-Pitaevskii equation

We consider Bose condensed $`F=1`$ spinor BEC’s with internal degrees of freedom for both ferromagnetic and antiferromagnetic cases. Here the order parameters are characterized by the hyperfine sublevels $`m_F=1,0,1`$. We start with the standard Hamiltonian by Ohmi and Machidaohmi , and Hoho :

$`\widehat{}_{rot}`$ $`=`$ $`\widehat{}{\displaystyle 𝑑𝐫𝛀\underset{j}{}\mathrm{\Psi }_j^{}(𝐫\times 𝐩)\mathrm{\Psi }_j}`$
$`\widehat{}`$ $`=`$ $`{\displaystyle }d𝐫[{\displaystyle \underset{ij}{}}\mathrm{\Psi }_i^{}\{h(𝐫)\mu _i\}\mathrm{\Psi }_j\delta _{ij}`$ (1)
$`+{\displaystyle \frac{g_\mathrm{n}}{2}}{\displaystyle \underset{ij}{}}\mathrm{\Psi }_i^{}\mathrm{\Psi }_j^{}\mathrm{\Psi }_j\mathrm{\Psi }_i`$
$`+{\displaystyle \frac{g_\mathrm{s}}{2}}{\displaystyle \underset{\alpha }{}}{\displaystyle \underset{ijkl}{}}\mathrm{\Psi }_i^{}\mathrm{\Psi }_j^{}(\widehat{F}_\alpha )_{ik}(\widehat{F}_\alpha )_{jl}\mathrm{\Psi }_k\mathrm{\Psi }_l].`$
Here

$`h(𝐫)={\displaystyle \frac{\mathrm{}^2^2}{2m}}+V(𝐫)`$ (2)
is one-body Hamiltonian. The quantity $`V(𝐫)=\frac{1}{2}m(2\pi \nu _r)^2(x^2+y^2)`$ is the external confinement potential such as an optical potential. The scattering lengths $`a_0`$ and $`a_2`$ characterize collisions between atoms through the total spin 0 and 2 channels, respectively, $`g_\mathrm{n}=\frac{4\pi \mathrm{}^2}{m}\frac{a_0+2a_2}{3}`$ is interaction strength through the “density” channel, and $`g_\mathrm{s}=\frac{4\pi \mathrm{}^2}{m}\frac{a_2a_0}{3}`$ is that through the “spin” channel. The subscripts $`\alpha =(x,y,z)`$ and $`i,j,k,l=(0,\pm 1)`$ correspond to the above three species. The chemical potentials for the three components $`\mu _i`$ $`(i=0,\pm 1)`$ satisfy $`\mu _1\mu _0=\mu _0\mu _1`$. We introduce $`\mu =\mu _0`$ and $`\mu ^{}=\mu _1\mu _0`$. The angular momentum operators $`\widehat{F_\alpha }(\alpha =x,y,z)`$ can be expressed in matrices as

$`F_x={\displaystyle \frac{1}{\sqrt{2}}}\left(\begin{array}{ccc}0& 1& 0\\ 1& 0& 1\\ 0& 1& 0\end{array}\right),`$ (6)
$`F_y={\displaystyle \frac{i}{\sqrt{2}}}\left(\begin{array}{ccc}0& 1& 0\\ 1& 0& 1\\ 0& 1& 0\end{array}\right),`$ (10)
$`F_z=\left(\begin{array}{ccc}1& 0& 0\\ 0& 0& 0\\ 0& 0& 1\end{array}\right).`$ (14)
Following the standard procedure, the extended Gross-Pitaevskii (GP) equation in rotation frame is obtained as

$`[\{h(𝐫)\mu _i`$ $`+`$ $`g_\mathrm{n}{\displaystyle \underset{l}{}}|\psi _l|^2\}\delta _{ij}`$ (15)
$`+g_\mathrm{s}{\displaystyle \underset{\alpha }{}}{\displaystyle \underset{kl}{}}\{(F_\alpha )_{ij}(F_\alpha )_{kl}\psi _k^{}\psi _l\}`$
$`i\mathrm{}𝛀\times 𝐫\delta _{ij}]\psi _j=0.`$
These coupled equations for the $`j`$-th condensate wave function $`\psi _j=\mathrm{\Psi }_j`$ $`(j=0,\pm 1)`$ are used to calculate various properties of vortices in the following. Here we take the external rotation as $`𝛀=\mathrm{\Omega }\widehat{𝒛}`$ and assume uniformity along z direction.

### II.2 Numerical procedure

The stationary states of the extended GP equation are defined as local minima of the energy functional

$`E[\psi _i,\psi _i^{}]={\displaystyle 𝑑𝐫\left[\underset{i}{}E_i(𝐫)+E_s(𝐫)\right]}\mu ^{}M𝛀𝐋,`$ (16)
where $`E_i`$ and $`E_s`$ are defined by

$`E_i(𝐫)`$ $`=`$ $`\psi _i^{}\left\{h(𝐫)\mu +{\displaystyle \frac{g_n}{2}}{\displaystyle \underset{k}{}}|\psi _k|^2\right\}\psi _i,`$ (17)
$`E_s(𝐫)`$ $`=`$ $`{\displaystyle \frac{g_s}{2}}{\displaystyle \underset{\alpha }{}}\left\{{\displaystyle \underset{k,l}{}}(\psi _k^{}(\widehat{F}_\alpha )_{kl}\psi _l)\right\}^2,`$ (18)
$`𝛀𝐋`$ denotes

$`𝛀𝐋=i\mathrm{}\mathrm{\Omega }{\displaystyle 𝑑𝐫\underset{i}{}\psi _i^{}\left(x\frac{}{y}y\frac{}{x}\right)\psi _i},`$ (19)
and $`M=𝑑𝐫_ii|\psi _i(𝐫)|^2`$ is the total magnetization.

The numerical algorithm used to minimize the energy functional in scalar BECcastin ; garcia ; aftalion ; feder  can be extended to the present system. Following this procedure, the initial $`\psi _j`$ given randomly are modified using the local gradient of the energy functional $`E[\psi _i,\psi _i^{}]`$ as

$`\psi _j(\tau +\mathrm{\Delta }\tau )=\psi _j(\tau ){\displaystyle \frac{\delta E[\psi _i,\psi _i^{}]}{\delta \psi _j^{}}}\mathrm{\Delta }\tau .`$ (20)
This equation means that the order parameters $`\psi _j`$, parameterized by a ‘fictitious time’ $`\tau `$, roll along the slope of the energy functional. Equation (20) is rewritten as

$`\mathrm{}_\tau \psi _j(\tau )`$ $`=`$ $`[\{h(𝐫)(\mu (\tau )+\mu ^{}j)+g_n{\displaystyle \underset{k}{}}|\psi _k|^2\}\delta _{jk}`$ (21)
$`+{\displaystyle \frac{g_s}{2}}{\displaystyle \underset{\alpha }{}}\left\{{\displaystyle \underset{l,p}{}}(\psi _l^{}(\widehat{F}_\alpha )_{lp}\psi _p)(\widehat{F}_\alpha )_{jk}\right\}`$
$`+i\mathrm{}\mathrm{\Omega }(x{\displaystyle \frac{}{y}}y{\displaystyle \frac{}{x}})\delta _{jk}]\psi _k(\tau ).`$
This is the Gross-Pitaevskii equation for imaginary times $`\tau =it`$. In each time step, $`\mu (\tau )`$ is adjusted to preserve the total number of particles in the system

$`N={\displaystyle \underset{j}{}}{\displaystyle 𝑑𝐫|\psi _j(𝐫)|^2}.`$ (22)
For $`\tau \mathrm{}`$, $`\psi _j`$ converges to the stationary state, corresponding to one of the local minima of the energy functional (16). For $`\tau =\mathrm{}`$ $`\psi _j`$ satisfies

$`{\displaystyle \frac{\delta E[\psi _i,\psi _i^{}]}{\delta \psi _j^{}}}|_\tau \mathrm{}=0,`$ (23)
and Eq. (21) becomes equivalent to the extended GP equation (15).

We take the initial state of each component as

$`\psi _j(𝐫,\tau =0)=\sqrt{n_{\mathrm{TF}}(𝐫)}\eta _j\mathrm{exp}[iS_j(𝐫)],`$ (24)
where $`n_{TF}(𝐫)`$ is the density profile within the Thomas-Fermi (TF) approximation:

$`n_{\mathrm{TF}}(𝐫)=\{\begin{array}{cc}\frac{\mu _{\mathrm{TF}}V}{g_n}\hfill & \text{for }g_s>0\\ \frac{\mu _{\mathrm{TF}}V}{g_n+g_s}\hfill & \text{for }g_s<0\end{array},`$ (27)
and $`\eta _j`$ represents the ratio of each component. The phase is given by

$`S_j(𝐫)={\displaystyle \underset{k}{}}w_j^{(k)}\theta _j^{(k)}+\alpha _j,`$ (28)
$`w_j^{(k)}`$ is the winding number of the $`j`$-th condensate, $`\theta _j^{(k)}`$ is the polar angle of the coordinate $`(x^{(k)},y^{(k)})`$ whose origin is located at the $`k`$-th vortex core, and $`\alpha _j`$ is relative phase between the three components.

It is convenient to describe the condensates in terms of the three components $`\psi _\alpha (\alpha =x,y,z)`$ where the quantization axis is taken along the $`\alpha `$ direction:

$`\left(\begin{array}{c}\psi _x(𝐫)\\ \psi _y(𝐫)\\ \psi _z(𝐫)\end{array}\right)={\displaystyle \frac{1}{\sqrt{2}}}\left(\begin{array}{ccc}1& 0& 1\\ i& 0& i\\ 0& \sqrt{2}& 0\end{array}\right)\left(\begin{array}{c}\psi _1(𝐫)\\ \psi _0(𝐫)\\ \psi _1(𝐫).\end{array}\right)`$ (38)
We then define a couple of real vectors as

$`𝒎=(m_x,m_y,m_z)=\mathrm{Re}(\psi _x,\psi _y,\psi _z),`$ (39)
$`𝒏=(n_x,n_y,n_z)=\mathrm{Im}(\psi _x,\psi _y,\psi _z).`$ (40)
The $`l`$-vector, which points the direction of the local magnetization, is defined as $`𝒍=𝒎\times 𝒏`$. The corresponding unit vector is denoted by $`\widehat{𝒍}`$.

### II.3 Calculated system

The actual calculations are carried out by discretizing the two-dimensional space into 51$`\times `$51 mesh. We have performed extensive search to find stable vortices, starting with various vortex configurations, covering a wide range of the ferromagnetic and the antiferromagnetic interaction strength, $`g_s/g_n=0.20.2`$, and examining various axisymmetric and non-axisymmetric vortices. See Ref.isoshima1 ; isoshima2  for the classification of possible vortices in the axisymmetric case. We use the following parameters: the mass of a 87Rb atom $`m`$=1.44$`\times 10^{25}`$kg, the trapping frequency $`\nu _r`$=200Hz, and the particle number per unit length along the $`z`$ axis $`n_z`$=2.0$`\times 10^3/\mu `$m. The results displayed here are for $`g_s/g_n=0.02`$ (ferromagnetic case) and $`g_s/g_n=0.02`$ (antiferromagnetic case). The external rotation frequency $`\mathrm{\Omega }`$ is normalized by the harmonic trap frequency.

## III vortex structure

The vortex configurations are characterized by the combination of the winding number $`w_j`$ of $`\psi _j`$ ($`j=0`$, $`\pm 1`$) denoted by $`w_1,w_0,w_1`$, where $`w_j`$ denotes the phase change by $`2\pi w_j`$ when the wave function goes around the phase singularity.

The spin term (18) of the total energy is rewritten as

$`E_s(𝐫)={\displaystyle \frac{g_s}{2}}\left[n^2(𝐫)\left|2\psi _1(𝐫)\psi _1(𝐫)\psi _0^2(𝐫)\right|^2\right],`$ (41)
where $`n(𝐫)=_j|\psi _j(𝐫)|^2`$ is the total density. The spin texture in the ground state without rotation is determined by this energy. By minimizing Eq.(41), the relative phases $`\alpha _j`$ are shown to satisfy

$`2\alpha _0=\alpha _1+\alpha _1+n\pi ,`$ (42)
where $`n`$ is an integer, and the odd (even) $`n`$ corresponds to the antiferromagnetic (ferromagnetic) situationisoshima1 .

### III.1 Axisymmetric vortex

It also follows from Eq. (41) that $`w_1,w_0,w_1`$ of axisymmetric vortices satisfies

$`2w_0=w_1+w_1.`$ (43)
Thus, possible candidates for the stable state are the non-vortex state $`0,0,0`$ and the vortex configurations: $`1,0,1`$, $`1,\times ,0`$, $`0,\times ,1`$, and $`1,1,1`$, which exhaust all the combinations of the winding numbers less than or equal to unity. The thermodynamic stability of these vortices are demonstrated in Ref.isoshima2 . Here we concentrate on the possibility of combinations with higher winding numbers, which will be shown to be stable only in the ferromagnetic case.

Figure 1 displays the density and $`l`$-vector profiles of a new continuous vortex $`0,1,2`$ found stable for the ferromagnetic interaction ($`g_s<0`$). It is seen that $`\psi _1`$ with zero winding number $`w_1=0`$ occupies the central region of the harmonic trap and $`\psi _1`$ with the higher winding number $`w_1=2`$ fill in the circumference region. The intermediate region is occupied by $`\psi _0`$ component which has a singularity $`w_0=1`$ at the center of the trap. The resulting total density is non-singular and have a smooth spatial variation described by a Gaussian form except for the outermost region. This vortex is equivalent to the topological structure called the Mermin-Ho vortex in superfluid 3Hemermin  and Skyrmion in generalstoof .

The axisymmetric $`0,1,2`$ vortex may be parametrized as

$`\left(\begin{array}{c}\psi _1(𝐫)\\ \psi _0(𝐫)\\ \psi _1(𝐫)\end{array}\right)=\sqrt{n(r)}\left(\begin{array}{c}\mathrm{cos}^2\frac{\beta }{2}\\ \sqrt{2}e^{i\varphi }\mathrm{sin}\frac{\beta }{2}\mathrm{cos}\frac{\beta }{2}\\ e^{2i\varphi }\mathrm{sin}^2\frac{\beta }{2}\end{array}\right)`$ (50)
where the bending angle $`\beta (r)`$ runs over $`0\beta (r)\pi `$ and $`\varphi `$ signifies the polar angle in polar coordinates. The spin direction is denoted by the $`l`$-vector and is given as $`𝒍(r)=\widehat{z}\mathrm{cos}\beta +\mathrm{sin}\beta (\mathrm{cos}\varphi \widehat{x}+\mathrm{sin}\varphi \widehat{y})`$ where $`\beta `$ varies from $`\beta (0)=0`$ to $`\beta (R)=\frac{\pi }{2}(=\pi )`$ for MH (Anderson-Toulouse (AT)) ($`R`$ is the outer boundary of the cloud). Thus the spin moment is flared out to the radial direction and at the circumference it points outward for MH and downwards for AT (for schematic $`l`$-vector structure, see Fig.18 in Ref.salomaa  ).

In Fig.2 we show the spatial dependence of the $`l_z`$-component along the radial direction, namely, the spatial dependence of the bending angle $`\beta (r)`$ for the MH vortex. As the magnetization $`M`$ decreases, the local magnetization in the condensate surface changes from positive to negative passing through zero. It means that the $`l`$-vector in this vortex flares out radially to orient almost horizontally $`\beta (r=R)=\frac{\pi }{2}`$ for $`M/N0.5`$ and to point downward for $`\beta (r=R)=\pi `$ for $`M/N0`$. The former (latter) corresponds literally to the Mermin-Ho (Anderson-Toulouse) vortex. This is simply because as $`M`$ decreases, the spin-down component $`\psi _1`$ with $`w=2`$ increases in the outer region. Thus we can control these MH and AT vortices by merely changing the total magnetization.

As pointed out in the previous papermizushima , however, the situation is completely different from the case of superfluid 3He-A where the stability of the MH vortex is due to the constraint that the $`l`$-vector be perpendicular to the vessel wallsalomaa . These vortex configurations in ferromagnetic BEC are created naturally under the condition of a given total number and magnetization, both of which are well controlled in a harmonic trap experiment.

In comparison with the $`0,1,2`$ vortex, the spin textures of other vortex configurations, such as the $`1,\times ,0`$ and $`1,0,1`$ vortices, have a different natureisoshima1 . In $`1,\times ,0`$ vortex, the spin moment is suddenly reversed near the vortex core of $`\psi _1`$ component because of the absence of $`\psi _0`$ component. In $`1,0,1`$ vortex it can vary continuously around the vortex core. In this configuration, however, since the condensate at the center of the trap consists only of the polar stateho , the spin texture has a singularity. Thus only $`0,1,2`$ vortex can have a non-singular and continuous spin texture under slow rotation.

It is easy to calculate the total angular momentum $`L_z`$ of the axisymmetric vortices; by using the total particle number and the total magnetization, it is simply written as

$`{\displaystyle \frac{L_z}{\mathrm{}N}}=w_0+w^{}{\displaystyle \frac{M}{N}},`$ (51)
where the total magnetization is written as $`M=𝑑𝐫n(𝐫)\widehat{l}_z`$ and we have introduced $`w_jw_0=jw^{}`$. Thus the spin textures with a net spin polarization carry the angular momentum, i.e. the superflow. For $`0,1,2`$ vortex, $`\frac{L_z}{\mathrm{}N}=1\frac{M}{N}`$. This simple formula has the following physical meaning. (i)At $`M=N`$, $`L_z=0`$ because $`\psi _1`$ has no winding. (ii)At $`\frac{M}{N}=\frac{1}{2}`$, $`\frac{L_z}{\mathrm{}N}`$ is exactly equal to $`\mathrm{}/2`$, corresponding to the MH vortex.

In the $`0,1,2`$ state with the higher winding, the non-winding component $`\psi _1`$ works as a “pinning potential” for the remaining $`\psi _0`$ and $`\psi _1`$, thereby making the state stable in the lower rotation drive. In particular, $`\psi _1`$ with $`w_1=2`$ is stabilized by the presence of the $`\psi _1`$ due to the ferromagnetic interaction. For a very small antiferromagnetic interaction ($`g_s=0.005g_n`$) and non-magnetic case ($`g_s/g_n=0`$), the vortex with the $`w_1=2`$ becomes unstable and splits into a couple of $`w_1=1`$ vortices (see Fig.3). This configuration is equivalent to the vortex found by Yip (phase IV in Ref.yip ) and is always unstable for the large $`g_s/g_n`$ ($`>0`$).

### III.2 Non-Axisymmetric vortex

To investigate the possibility of non-axisymmetric vortices, let us first recapitulate the axisymmetric $`1,1,1`$ vortex. In the axisymmetric case, the total density of the $`1,1,1`$ vortex is equivalent to the one in scalar BEC, having the singularity at the center of the trap where the potential energy is minimum. The axisymmetric singular $`1,1,1`$ vortex is always unstable even in the higher rotationisoshima2 . However, by displacing the vortex cores of each component from the center of the trap, the $`1,1,1`$ vortex can be stable as a non-axisymmetric non-singular type in the lower rotation frequency. This works favorably to gain more condensation energy at the center of the trap. A similar situation is seen near the $`\mathrm{\Omega }_{c2}`$ in the homogeneous system, where the singular vortex lattice, called the Abrikosov lattices, is never favored in the entire antiferromagnetic regionkita .

For the antiferromagnetic case, $`\psi _1`$ and $`\psi _1`$ overlap to minimize the spin-dependent energy $`E_s(𝐫)`$. In Fig.4, we show the density profiles and $`g_s`$-dependence of two different $`1,1,1`$ vortices. In the non-magnetic limit ($`g_s/g_n=0`$), the two states are completely equivalent. As the spin interaction $`g_s/g_n`$ increases, striking differences grow between the two states. For $`g_s/g_n0.4`$, the vortex cores of the state presented in the left of Fig.4 (b) (the split-(I) state) collapse and the amplitude of each order parameter cannot recover near the vortex cores, i.e. this state behaves like the phase separation in $`x`$-$`y`$ plane. In contrast, the state shown in the right of Fig.4 (b) (the split-(II) state) forms the regular cores. As a result, the split-(II) state is energetically favorable over the split-(I) state for the antiferromagnetic situation. In very small spin interaction range ($`g_s0.02g_n`$), which is a realistic parameter, however, the split-(I) state can be favored over the split-(II), though the energies of two states are very close to each other.

A third vortex configuration is displayed in Fig.5 where vortex cores of each component are displaced to form a triangle. In the non-magnetic limit, the three components are completely equivalent and the three singularities form a regular triangle. As the spin interaction $`g_s`$ increases, the $`\psi _1`$ and $`\psi _1`$ components overlap locally because of the antiferromagnetic interaction and this vortex configuration starts to deform from a regular triangle shape. For larger spin interactions ($`g_s>0.4`$), this vortex becomes unstable due to the overlap between $`\psi _1`$ and $`\psi _1`$ components.

Figure 6 (a) shows the density profiles of a non-axisymmetric non-singular $`1,1,1`$ vortex for the ferromagnetic case. Two singularities of $`\psi _1`$ and $`\psi _1`$ are displaced from the center. The $`\psi _0`$ component with the singularity at the center of the trap plays the role to prevent the phase separation favored in the ferromagnetic spin interaction. The spin texture in this state is displayed in Figs.6 (b) and (c) where the spin moments flip at the center of the trap.

In comparison with axisymmetric types, these non-axisymmetric vortices have the advantage that they can easily adjust themselves for a change in $`\mathrm{\Omega }`$. As $`\mathrm{\Omega }`$ increases, the two or three separate singularities adjust their mutual distance from the center and change the value of $`L_z`$ to gain the energy from the term $`\mathrm{\Omega }L_z`$. In this sense, the non-axisymmetric vortices are flexible against a change in $`\mathrm{\Omega }`$ compared with the axisymmetric ones.

## IV phase diagram

The phase diagrams are shown in a plane of the external rotation and the total magnetization by comparing the energies of various vortex configurations:

$`E={\displaystyle 𝑑𝐫\left[\underset{i}{}E_i(𝐫)+E_s(𝐫)\right]}\mathrm{\Omega }L_z.`$ (52)
It is noted that the region $`\mathrm{\Omega }<0.38`$ considered here corresponds to the single vortex region in the scalar BEC caseisoshima3 . Thus, the present single-vortex consideration may also be justified for $`\mathrm{\Omega }<0.38`$.

The resulting phase diagram is displayed in Fig.7 for (a) the ferromagnetic case and (b) the antiferromagnetic case. For the ferromagnetic case, a large area of the $`\mathrm{\Omega }M`$ plane is occupied by the $`0,1,2`$ vortex, including MH and AT. The non-axisymmetric $`1,1,1`$ vortex and the $`1,0,1`$ vortex are stabilized near $`M=0`$ and $`M=N`$, respectively. We find a large empty area in the intermediate $`M/N`$ region where neither single-vortex nor vortex-free states are stabilized at all because the phase separation in the ferromagnetic case prevents forming a uniform mixture of the three components even when the circulation is absent in the vortex-free state.

In the phase diagram for the antiferromagnetic case, in contrast, everywhere is occupied by a stable phase. This result is consistent with Fig.2(a) of Ref.isoshima2  over a wide range, except for the presence of two non-axisymmetric types near $`M0`$, i.e. $`1,1,1`$ split-(I) and triangle vortex. The $`1,1,1`$ triangle vortex is energetically indistinguishable from the phase-I vortex obtained by Yipyip . The phase-IV vortex given in Fig,3 of Ref.yip  is found unstable for $`g_s=0.02g_n`$ and does not appear in this phase diagram.

It is found for both cases that the stabilities of the non-axisymmetric vortices are restricted in a narrow region. On the other hand, the axisymmetric vortices are stable in a large area. As discussed in Section III, since finite density of a component at the vortex cores of the others supports its stability, the non-axisymmetric vortex becomes unstable in increasing $`M`$ where the number of $`\psi _1`$ grows while the others shrink.

Figure 8 shows the $`\mathrm{\Omega }`$-dependence of the angular momentum at $`M=0`$ for the two cases. For axisymmetric types, as shown in Eq.(51), the angular momentum of the system $`L_z`$ is fixed for a given $`M`$. Thus there is the need of changing the winding combinations $`w_1,w_0,w_1`$ so as to increase $`L_z`$, which means that the axisymmetric vortices do not have the adaptability for changes in $`\mathrm{\Omega }`$.

## V Discussion

The phase-separated state with $`w=0`$ is expected to be stable in a large empty region of Fig.7 (a). In this state, $`\psi _1`$ and $`\psi _1`$ components phase-separate along $`z`$-direction due to the ferromagnetic interaction. Namely, an arbitrary $`x`$-$`y`$ cross-section consists of only $`\psi _1`$ or $`\psi _1`$ component, and the spin-polarized state with $`M/N=\pm 1`$ is piled up along the z-direction. Neglecting the contribution from the boundary layer associated with the phase separation, we can estimate that the energy of this phase-separated state with $`w=0`$ and $`w=1`$ is equal to the energy of the scalar BEC with $`w=0`$ and $`w=1`$, respectively.

We compare in Fig.9 (a) the free energies of the $`0,1,2`$ vortex state and the phase-separated state. As shown in Fig.9 (a), the $`0,1,2`$ vortex with the three components is energetically favored over the phase-separated state, where the energy of the phase-separated state is given by the energy of the $`0,1,2`$ vortex at $`M/N=1`$, i.e. the scalar BEC with $`w=0`$. Thus, the composite state of the three components, such as the MH vortex, becomes “locally” stable under a rotation drive while the composite state may have “global” stability. It is possible to perform a similar discussion for the phase-separated state in higher rotation. The result of Fig.9 (b) shows that the $`1,0,1`$ vortex is favored over the phase-separated state with the winding, whose total density corresponds to the conventional singular vortex in the scalar BEC.

## VI Conclusion

We have presented possible vortex structures and the vortex phase diagram in the plane of external rotation $`\mathrm{\Omega }`$ and the total magnetization $`M/N`$ for the both cases of antiferromagnetic ($`g_s=0.02g_n`$) and ferromagnetic ($`g_s=0.02g_n`$) interaction. We have investigated the thermodynamic stability of the possible vortex configurations by solving the extended Gross-Pitaevskii equation for the spinor BEC with $`F=1`$.

For the ferromagnetic case, the stability of the continuous vortex, called the Mermin-Ho and Anderson-Toulouse vortex, is demonstrated, but these topological structures are found never stable under no rotation drive. Furthermore, these vortices can exist in the intermediate process (see Fig.5 in the Ref.isoshima4 ) proposed by Isoshima et al.isoshima4 ; nakahara , i.e. it may be created by making use of spin texture.

We have also discovered a couple of new non-axisymmetric vortices besides the two vortices found by Yipyip  for the antiferromagnetic case. The conventional singular vortex is found to be never favored in spinor BEC for both cases. It means that the total density profile is always non-singular and continuous. Therefore, the experimental procedure to image the magnetic patterns for each vortex configuration will be required as a special technique to identify these vortices.

###### Acknowledgements.

The authors thank T. Ohmi and T. Isoshima for useful discussions. One of the authors (TM) would like to acknowledges the financial support of Japan Society for the Promotion of Science (JSPS) for Young Scientists.