format: text

text:
# Normal modes and stability of phase separated trapped Bose-Einstein condensates

## I Introduction

There is much recent interest in Bose-Einstein condensation (BEC) in trapped gases as well as quantum phase transitions (QPT). Examples of the QPT include the Wigner electron solid melting transition, the Mott-Hubbard metal-insulator transition, and different magnetic transitions. In physical phenomena involving BECs quantum mechanics play a crucial role. In this paper we propose a new kind of quantum phase transition in phase-separated mixtures of BECs. In this transition, the distribution of the two components changes from a symmetric to an asymmetric shape. To explore the transition, we first investigate the stability of the symmetric phase by studying its normal modes. We find interface modes that become soft. These unstable modes can be caused by changing the trapping potential and are analogous to gravitational waves at the interface between two immiscible fluids. When the lowest frequency becomes zero, the instability sets in which determines the stability limit of the symmetric phase. We determine the actual phase boundary by comparing the energy between the symmetric phase and the asymmetric phase and find that the actual phase boundary and the instability boundary is not the same. This suggests that the transition is first order. The system may be a good laboratory to study issues of quantum metastability and tunnelling. The symmetric to asymmetric transition is the simplest quantum phase transition that one can imagine. Careful study of this problem should provide us new insight into this burgeoning field of discovery.

Systems of multicomponent condensates were first realized by the JILA group in a magnetic trap in 87Rb  and subsequently in an optical trap in 23Na (spinor condensate) . Experimentally at low fields, the spin exchange process can occur in an optically trapped condensate, leading to spin domains  with metastable behavior . Binary condensates in two hyperfine levels of 87Rb have been studied , most notably realizing a system of interpenetrating Bose fluids , measurements of phase dispersion , and a vortex state in a dilute-gas BEC . Theoretical treatment of such systems began in the context of superfluid helium mixtures  and spin-polarized hydrogen , and now has been extended to BEC in trapped gases .

The equilibrium density distributions of segregated mixtures in the absence of gravity have been studied numerically for different system parameters. Two types of configurations have been discussed: a symmetric  configuration, for which one component is inside the other one, and an asymmetric one in which the two components occupy the left and the right hand side of a sphere . Several theoretical papers have described collective excitations of multicomponent condensates. Busch et al.  studied collective excitations in the limit of weak interactions based on a trial function approach. Collective excitations of binary mixtures in the Thomas-Fermi (TF) limit have been investigated by Graham and Walls . Esry and Greene  based on the Hartree-Fock and random-phase approximation have numerically calculated the low-lying excitations of double BECs in a time-averaged orbiting potential trap in which gravity separates the centers of the two components. Gordon and Savage  have studied excitation spectrum as a function of number of particles. Mazets  has performed an analysis of wave dispersion on a boundary between two weakly segregated untrapped BECs.

Under a change of the interaction strength the two component BECs can undergo different types of macroscopic quantum phase transitions. For example, an increase in interspecies interactions $`a_{12}`$ results in a transition from a binary mixture to a phase separated state. Such phase separated condensates are analogous to a system of immiscible fluids with positive surface tension. Pu and Bigelow  numerically studied the frequency of collective excitations as a function of interspecies interactions $`a_{12}`$ and found appearance of imaginary normal mode frequencies at the phase separation transition point. Such imaginary modes describe the nature of system’s instability. The trapping potential is not crucial for this type of phase transition and similar imaginary modes also exist in the homogeneous (untrapped) BECs .

However, in the phase separated regime the trapping potential changes the system’s symmetry and results in a new type of quantum phase transition. Such transition occurs in trapped condensates upon the change of relative intraspecies interaction strength $`a_{11}/a_{22}`$. When this ratio differs significantly from one the less repulsive component is in the middle of the trap and the more repulsive component comprises outer shell. One can make the initial configuration unstable by manipulating the ratio $`a_{11}/a_{22}`$ close to one, e.g., by means of Feschbach resonances or by changing the ratio of the trapping frequencies. In this paper we study the collective excitations of the trapped phase separated condensates as a function of $`a_{11}/a_{22}`$ and show that normal modes with imaginary frequencies appear when this ratio approaches one. Effect of the interface tension results in a shift of the transition point of the symmetric-asymmetric transition in the two component BEC. In general, the asymmetric phase possesses a lower interface energy. On the other hand, since the degree of self-repulsion may differ between the two species, the less self-repulsive component will prefer to remain where the density is higher, while the other component moves to the low density regions outside. This favors the symmetric phase. Depending on the system parameters, one of these two energetic considerations will win out. These system parameters can be adjusted by changing the trapping frequencies, the relative particle numbers of the two species, and the interaction between the particles with Feschbach resonances. We first address the stability of the symmetric phase.

## II Basic formalism

Let us consider the two-component BEC in a spherically symmetric trap. The dynamics of the system is described by time dependent Gross-Pitaevskii equations

$$i\mathrm{}\frac{\mathrm{\Psi }_1}{t}=\frac{\mathrm{}^2}{2m}\mathrm{\Delta }\mathrm{\Psi }_1+V_{\mathrm{tr}}\mathrm{\Psi }_1+\frac{4\pi \mathrm{}^2}{m}(a_{11}|\mathrm{\Psi }_1|^2+a_{12}|\mathrm{\Psi }_2|^2)\mathrm{\Psi }_1,$$
(1)
$$i\mathrm{}\frac{\mathrm{\Psi }_2}{t}=\frac{\mathrm{}^2}{2m}\mathrm{\Delta }\mathrm{\Psi }_2+V_{\mathrm{tr}}\mathrm{\Psi }_2+\frac{4\pi \mathrm{}^2}{m}(a_{22}|\mathrm{\Psi }_2|^2+a_{12}|\mathrm{\Psi }_1|^2)\mathrm{\Psi }_2,$$
(2)
where $`\mathrm{\Psi }_{1,2}`$ are the condensate wave functions, $`V_{\mathrm{tr}}=m\omega _0^2r^2/2`$ is the trapping potential, $`\omega _0`$ is the trapping frequency, $`r`$ is the radial spherical coordinate, $`a_{ij}>0`$ are $`s`$wave scattering lengths. Here we suppose that trapping frequencies and particle masses are the same for both species. We discuss the general case in the last section.

We shall assume in this paper that the condition

$$a_{12}^2a_{11}a_{22}>0$$
(3)
is satisfied and, therefore, the condensates are phase-segregated. The first demonstration of a condensate binary mixture by Myatt et. al.  produced overlapping condensates of the $`|F=1,m_f=1>`$ and $`|F=2,m_f=2>`$ spin states of 87Rb. These states, however, possess different magnetic moments and, hence, the condensates experience different potentials in a magnetic trap which results in unequal displacement from the trap center by gravity. Later JILA experiments were performed on mixtures of $`|1,1>`$ and $`|2,1>`$ states . These two states have essentially identical magnetic moments, and feel identical confining potentials. In 87Rb, the scattering lengths for $`|1,1>`$ and $`|2,1>`$ states are known to be in the proportion $`a_{11}:a_{12}:a_{22}::1.03:1:0.97`$, with the average of the three being $`55(3)`$Å . Hence, the condition (3) is valid for the JILA experiments and the condensates are in a weakly segregated phase . One should mention that for earlier experiments with the $`|1,1>`$ and $`|2,2>`$ rubidium states the scattering lengths are $`a_{11}:a_{12}:a_{22}::1.007:1:1.01`$ . For such states the condition (3) is not fulfilled and the mixture exhibits behavior of a miscible system.

We suppose the condensates can be well described by the TF approximation. In this regime, the phase-segregated condensates overlap over the length scale $`\mathrm{\Lambda }=\xi /\sqrt{a_{12}/\sqrt{a_{11}a_{22}}1}`$, where $`\xi `$ is the healing length . For the JILA experiments on phase-segregated states $`\mathrm{\Lambda }47\xi `$. If, however, the penetration depth $`\mathrm{\Lambda }R`$, where $`R`$ is the size of the system, the condensates can be approximately treated as nonoverlapping, which we assume to be the case. The effect of the overlapping results in finite surface tension and can be included via boundary conditions at the interface. However, in the TF limit such effect is small, we shall discuss it in Sec. IV. If the condensates do not overlap one can neglect the last terms in Eqs. (1), (2). As a result, the dynamical equations for $`\mathrm{\Psi }_1`$ and $`\mathrm{\Psi }_2`$ decouple. However, the two condensate components are coupled by boundary conditions at the interface which require continuity of pressure and the normal velocity.

For definiteness we assume that stationary configuration is spherically symmetric with the central core dominated by the first component and an outer shell from the second species (see Fig. 1). The stationary density distribution $`n_i=|\mathrm{\Psi }_i|^2`$ of two components is given by

$$n_1=\frac{\mu _1}{G_{11}}\left(1\frac{r^2}{R_1^2}\right),0<r<R_{},$$
(4)
$$n_2=\frac{\mu _2}{G_{22}}\left(1\frac{r^2}{R_2^2}\right),R_{}<r<R_2,$$
(5)
where $`G_{ii}=4\pi \mathrm{}^2a_{ii}/m`$, $`R_i=\sqrt{2\mu _i/m\omega _0^2}`$. The normalization condition $`n_i𝑑V=N_i`$, where $`N_i`$ are the numbers of condensate particles, determines the chemical potentials $`\mu _i`$. The position of the phase boundary $`R_{}`$ is given by the condition that pressures exerted by both condensates are equal :

$$R_{}=R_2\sqrt{\frac{1\kappa \lambda }{(1\kappa )\lambda }},$$
(6)
where we have introduced dimensionless parameters

$$\kappa =\sqrt{\frac{a_{11}}{a_{22}}},\lambda =\frac{\mu _2}{\mu _1}.$$
(7)
Fig. 1. Symmetric configuration of BEC components.

The symmetric configuration is favorable when $`\kappa `$ differs from unity, with the less repulsive component being in the middle ($`a_{11}<a_{22}`$, that is $`\kappa <1`$). At the interface between the two components the stationary condensate densities are

$$n_1=\frac{\kappa (\mu _2\mu _1)}{G_{11}(1\kappa )},n_2=\frac{(\mu _2\mu _1)}{G_{22}(1\kappa )}.$$
(8)
$`n_2/n_1=\kappa <1`$.

One can rewrite the time dependent Gross-Pitaevskii equations in a hydrodynamic form which shows an analogy between our problem and the motion of two immiscible fluids. In the strong phase-segregated regime the dynamics of each components is described by the following hydrodynamic equations

$$\frac{n}{t}+(n𝐕)=0,$$
(9)
$$\frac{1}{2}mV^2+V_{\mathrm{tr}}\frac{\mathrm{}^2}{2m}\frac{1}{\sqrt{n}}\mathrm{\Delta }\sqrt{n}+Gn+m\frac{\mathrm{\Phi }}{t}=\mu ,$$
(10)
where $`𝐕`$ is the condensate velocity and $`\mathrm{\Phi }`$ is the velocity potential, $`𝐕=\mathrm{\Phi }`$. Eq. (9) is continuity equation for compressible flow, while Eq. (10) is the analog of the Bernoulli equation with $`(e+P)/n=\frac{\mathrm{}^2}{2m}\frac{1}{\sqrt{n}}\mathrm{\Delta }\sqrt{n}+Gn`$, where $`e`$ is the internal energy density and $`P`$ is the pressure. The kinetic energy (quantum) pressure is omitted in the TF limit and, therefore, $`e=Gn^2/2`$. As a result, from Eq. (10) we obtain for the pressure

$$2P=Gn^2=\mu nnV_{\mathrm{tr}}mn\frac{\mathrm{\Phi }}{t}\frac{1}{2}mnV^2.$$
(11)
The trapping potential plays the role of gravitational potential in hydrodynamics. The linearized hydrodynamic equations determine the condensate normal modes. In terms of perturbation in the condensate density $`n^{}`$ and the velocity potential $`\mathrm{\Phi }`$ the normal-mode amplitudes satisfy the coupled equations

$$i\omega n^{}=(𝐕n^{})+(n\mathrm{\Phi }),$$
(12)
$$i\omega \mathrm{\Phi }=𝐕\mathrm{\Phi }+\frac{G}{m}n^{}\frac{\mathrm{}^2}{4m^2n}\left[n\left(\frac{n^{}}{n}\right)\right],$$
(13)
where $`n`$ is the static condensate density and $`𝐕`$ is the static condensate velocity. In our problem $`𝐕=\mathrm{𝟎}`$, hence, $`n^{}i\omega m\mathrm{\Phi }/G`$ and the equation for the perturbation in velocity potential is

$$\frac{2\omega ^2}{\omega _0^2}\mathrm{\Phi }2r\frac{\mathrm{\Phi }}{r}+(R^2r^2)\mathrm{\Delta }\mathrm{\Phi }=0.$$
(14)
One can seek solution of this equation in the form

$$\mathrm{\Phi }=\mathrm{\Phi }(r)Y_{lm}(\theta ,\varphi ),$$

where $`Y_{lm}(\theta ,\varphi )`$ is the spherical harmonics. Then, using $`\mathrm{\Delta }Y_{lm}(\theta ,\varphi )=l(l+1)/r^2`$, we obtain the following equation for the radial function $`\mathrm{\Phi }(r)`$

$$r^2(r^2R^2)\frac{d^2\mathrm{\Phi }}{dr^2}+(4r^22R^2)r\frac{d\mathrm{\Phi }}{dr}+\left[l(l+1)R^2r^2\left(\frac{2\omega ^2}{\omega _0^2}+l(l+1)\right)\right]\mathrm{\Phi }=0.$$
(15)
Let us derive boundary conditions for $`\mathrm{\Phi }`$ at the interface. The change in density at the interface is the sum of the term from the movement of the interface and that from the perturbation of the density. Let $`\varsigma =\varsigma (\theta ,\varphi ,t)`$ be a small deviation of the $`r`$coordinate of the interphase from its stationary value $`R_{}`$. Then the condensate density at the interface is $`n(R_{}+\varsigma ,t)`$ $`n(R_{},t)+`$ $`\varsigma n(R_{},t)/r`$, and, hence, the pressure is given by

$$P=\frac{Gn^2}{2}\frac{Gn^2(R_{},t)}{2}+Gn(R_{},t)\frac{n(R_{},t)}{r}\varsigma .$$
(16)
In terms of the density perturbation $`n^{}`$, $`n(R_{},t)=n(R_{})+n^{}(R_{},t)`$, we obtain

$$P\frac{Gn^2(R_{})}{2}+Gn(R_{})n^{}(R_{},t)n(R_{})m\omega _0^2R_{}\varsigma .$$
(17)
Continuity of the pressure at the interface results in the following boundary condition

$$G_{11}n_1(R_{})n_1^{}(R_{},t)n_1(R_{})m\omega _0^2R_{}\varsigma _1=G_{22}n_2(R_{})n_2^{}(R_{},t)n_2(R_{})m\omega _0^2R_{}\varsigma _2.$$

Taking into account $`Gn^{}=im\omega \mathrm{\Phi }`$ and $`\varsigma /t=V_r=\mathrm{\Phi }/r`$ we finally obtain the following boundary condition for the velocity potential at the interface

$$\omega _0^2R_{}\frac{\mathrm{\Phi }_1}{r}\omega ^2\mathrm{\Phi }_1=\kappa \omega _0^2R_{}\frac{\mathrm{\Phi }_2}{r}\kappa \omega ^2\mathrm{\Phi }_2.$$
(18)
Continuity of $`V_r`$ gives another boundary condition

$$\frac{\mathrm{\Phi }_1}{r}=\frac{\mathrm{\Phi }_2}{r}.$$
(19)
Eq. (15) and the boundary conditions (18), (19) compose a complete set of equations necessary to determine normal modes of the system.

## III Normal modes

For one component condensates the solutions of Eq. (15) which are finite at $`r=0`$ and $`r=R`$ are

$$\mathrm{\Phi }(r)r^lP_n^{(l+\frac{1}{2},0)}\left(1\frac{2r^2}{R^2}\right),$$
(20)
where $`P_n^{(l+\frac{1}{2},0)}(x)`$ are Jacobi polynomials, $`n`$ is the radial quantum number. The corresponding eigenfrequencies are given by

$$\omega ^2=\omega _0^2(l+n(2n+2l+3)).$$
(21)
For the two component condensate the normal mode frequencies are different from the general case and new modes appear. However, we shall show that frequencies corresponding to $`n=0`$ remain the same as for the one component system. Rewriting the function $`\mathrm{\Phi }=r^l\eta (r^2/R^2)`$ in Eq. (15) results in a hypergeometric equation for $`\eta (\xi )`$

$$\xi (\xi 1)\eta _{\xi \xi }^{\prime \prime }+\left[\left(l+\frac{5}{2}\right)\xi l\frac{3}{2}\right]\eta _\xi ^{}+\frac{1}{2}\left(l\frac{\omega ^2}{\omega _0^2}\right)\eta =0.$$
(22)
Hence, general solutions for the inner and outer condensates which are regular at $`r=0`$ and $`r=R_2`$ respectively have the form

$$\mathrm{\Phi }_1=C_1r^lF(\alpha ,\beta ,l+3/2,\frac{r^2}{R_1^2}),\mathrm{\Phi }_2=C_2r^lF(\alpha ,\beta ,1,1\frac{r^2}{R_2^2}),$$
(23)
where $`F`$ is hypergeometric function, $`C_1`$, $`C_2`$ are constants and

$$\alpha =\frac{1}{2}\left[l+\frac{3}{2}\sqrt{l^2+l+\frac{9}{4}+\frac{2\omega ^2}{\omega _0^2}}\right],\beta =\frac{1}{2}\left[l+\frac{3}{2}+\sqrt{l^2+l+\frac{9}{4}+\frac{2\omega ^2}{\omega _0^2}}\right].$$

Using Eq. (23), the boundary conditions (18), (19) and the mathematical identity

$$\frac{d}{dz}F(\alpha ,\beta ,\gamma ,z)=\frac{\alpha \beta }{\gamma }F(\alpha +1,\beta +1,\gamma +1,z),$$

we obtain the equation for the normal mode frequencies

$$\frac{\omega ^2}{\omega _0^2}=(1\kappa )\frac{\left[l(l+3/2)+(l\omega ^2/\omega _0^2)\lambda xs_1(\omega ,x)\right]\left[(l\omega ^2/\omega _0^2)xs_2(\omega ,x)l\right]}{\left[l(l+3/2)(\kappa 1)+x(l\omega ^2/\omega _0^2)\left(\kappa \lambda s_1(\omega ,x)+(l+3/2)s_2(\omega ,x)\right)\right]},$$
(24)
where

$$s_1(\omega ,x)=\frac{F(\alpha +1,\beta +1,l+5/2,\lambda x)}{F(\alpha ,\beta ,l+3/2,\lambda x)},s_2(\omega ,x)=\frac{F(\alpha +1,\beta +1,2,1x)}{F(\alpha ,\beta ,1,1x)},x=\frac{R_{}^2}{R_2^2}.$$

One of the solutions of Eq. (24) is $`\omega ^2=l\omega _0^2`$, which coincide with those for one component condensate. For this solution the components oscillate in-phase and $`\mathrm{\Phi }_1=\kappa \mathrm{\Phi }_2r^lY_{lm}(\theta ,\varphi )`$. Another exact solution is $`\omega ^2=5\omega _0^2`$ with $`\mathrm{\Phi }_115r^2/3R_1^2`$, $`\mathrm{\Phi }_2\lambda (15r^2/3R_2^2)`$ which corresponds to $`l=0`$, $`n=1`$. For this solution the components oscillate out-of-phase if $`\sqrt{3/5}R_1<R_{}<\sqrt{3/5}R_2`$ and in-phase otherwise.

### A Numerical results

We solve Eq. (24) numerically and find normal mode frequencies $`\omega `$ of the two component condensate as a function of the parameter $`\kappa =\sqrt{a_{11}/a_{22}}`$ for different fixed ratios $`R_{}/R_2`$. The ratio $`R_{}/R_2`$ can be directly measured experimentally. Fig. 2 shows the normal mode frequencies for $`l=0`$. At fixed $`l`$ there is an infinite number of branches which correspond to different radial quantum numbers $`n`$. We plot only the lowest modes with $`n=0,1,2`$. In Fig. 3 the normal modes are estimated for $`l=1`$. In the limiting case $`a_{11}=a_{22}`$ ($`\kappa =1`$) the condensates behave as an one component system and the condensate modes coincide with Stringari’s result (21). Fig. 4 shows the low frequency modes that become imaginary at $`\kappa >1`$. These modes are peculiar for two component systems and are analogous to the waves at the interface between two layers of immiscible fluids under gravity . In trapped condensates the gradient of trap potential plays the role of gravitational field. As soon as $`\kappa `$ becomes greater than $`1`$, which means $`n_2>n_1`$ at the interface, the system becomes unstable the same way as two immiscible fluids in gravitational field when the more dense layer is on the top. One should mention that normal modes with imaginary frequency also appear at the transition point from a binary mixture into a phase separated state which occurs at $`a_{12}=\sqrt{a_{11}a_{22}}`$ . In that case, the fastest decaying mode occurs at a finite wavevector, resulting in a quantum spinodal transition .

Fig. 2. Normal mode frequencies for $`l=0`$ and $`n=1,2`$ as a function of $`\kappa =\sqrt{a_{11}/a_{22}}`$. The mode with $`n=2`$ depends on position of the interface, the mode frequency is estimated for $`R_{}/R_2=0.25,0.5,0.75`$. For $`\kappa =1`$ the mode frequencies coincide with those for the one component condensate.

Fig. 3. Normal mode frequencies for $`l=1`$ and $`n=0,1,2`$.

Fig. 4. Low frequency modes as a function of $`\kappa =\sqrt{a_{11}/a_{22}}`$ for different $`l=1,2,3`$. The position of the interface is $`R_{}=R_2/2`$. Frequencies become imaginary at $`\kappa >1`$.

### B Low frequency modes

Here we study the low frequency modes in detail. In the region $`|1\kappa |1`$the mode frequencies are small: $`|\omega |\omega _0`$. In this limit one can take $`\omega =0`$ in the right side of Eq. (24) and put $`\kappa ,\lambda 1`$ in the multiple. As a result, we obtain

$$\omega ^2\omega _0^2(1\kappa )f(l,x),$$
(25)
where

$$f(l,x)=\frac{l\left(l+3/2+xs_1(0,x)\right)\left(xs_2(0,x)1\right)}{x\left[s_1(0,x)+(l+3/2)s_2(0,x)\right]}.$$

Eq. (25) describes behavior of the low frequency modes in the region close to the point of instability $`\kappa =1`$. In this region $`\omega \sqrt{1\kappa }`$ and becomes imaginary when $`\kappa >1`$. In Fig. 5 we plot $`f`$ as a function of $`x=R_{}/R_2`$ for different $`l`$. Imaginary part of frequencies is greater for larger $`l`$ and decreases with increasing $`R_{}/R_2`$, that is bigger inner droplets are more stable. Figs. 6, 7 show the radial distribution of the density perturbation for the low frequency modes with different $`l`$. In our estimates we choose the radius of the inner droplet $`R_{}=R_2/2`$. The density perturbations are normalized so that $`n_2^{}(R_2)=1`$. For the low frequency modes the two condensates oscillate out-of-phase. Fig. 6 corresponds to $`\kappa =0.9`$ (stable region), while for Fig. 7 $`\kappa =1.1`$ (unstable regime). One can see that the normal mode profile undergoes no changes at the point of phase transition, although the mode frequencies become imaginary at this point. However, due to the relation $`n^{}i\omega m\mathrm{\Phi }/G`$, the phase shift between oscillations of $`n^{}`$ and $`\mathrm{\Phi }`$ changes at the instability point from $`\pi /2`$ to $`0`$. Modes with small $`l`$ are delocalized and the whole condensate is involved in oscillations. With increasing $`l`$, however, the modes become localized near the interface and in the limit $`l\mathrm{}`$ they are similar to gravitational waves at the surface of deep water with the dispersion relation $`\omega \sqrt{k}`$.

Fig. 5. Dependence of the factor $`f`$ in Eq. (25) on the position of the interface $`R_{}/R_2`$ for modes with different $`l=1,2,3`$.

Fig. 6. Perturbation in the condensate density of the two components $`n_1^{}(r)`$ and $`n_2^{}(r)`$ for the low frequency modes with different $`l=1,2,3`$. The position of the interface is $`R_{}=R_2/2`$ (dash line) and $`\kappa =0.9`$ (stable region).

Fig. 7. The same as in Fig. 6, but for $`\kappa =1.1`$ (unstable region).

## IV Effect of interface tension

Interface tension is small in the TF limit and results in corrections of the order of $`\xi /R_2`$ to the normal mode frequencies. However, for the low frequency modes in the vicinity $`\kappa 1`$ the effect of the interface tension is substantial because the mode frequencies themselves are close to zero. The surface tension $`\sigma `$ modifies the boundary condition for the pressure so that the pressure difference at the interface is equal to the surface tension pressure $`P_1P_2=\sigma (1/r_1+1/r_2)`$, where $`r_1`$, $`r_2`$ are the principle radii of curvature. The stationary position of the interface is determined by the equation

$$\frac{G_{11}n_1^2}{2}=\frac{G_{22}n_2^2}{2}+\frac{2\sigma }{R_{}},$$

which results in the following ratio of densities at the interface

$$\frac{n_2}{n_1}=\sqrt{\kappa ^2\frac{4\sigma }{G_{22}n_1^2R_{}}}=\kappa _{\mathrm{eff}}.$$

The interface shape oscillates during the condensate motion which produces oscillations of the surface pressure. For interface shape close to spherical

$$\frac{1}{r_1}+\frac{1}{r_2}=\frac{2}{R_{}}\frac{2\varsigma }{R_{}^2}\frac{\mathrm{\Delta }_{\theta ,\varphi }\varsigma }{R_{}^2},$$
(26)
where $`\mathrm{\Delta }_{\theta ,\varphi }`$ is the angular part of Laplace’s operator. Also one should take into account that $`\sigma `$ itself depends on local density and varies during the interface motion. We assume $`\sigma \sqrt{\frac{1}{G_{11}n_1}+\frac{1}{G_{22}n_2}}n_1n_2`$ . Then, taking into account $`n(R_{}+\varsigma ,t)`$ $`n(R_{})+n^{}(R_{},t)+`$ $`\varsigma n(R_{})/r`$, we obtain

$$\sigma =\sigma _0\left(1+\frac{m}{G_{11}n_1}\left(i\omega \mathrm{\Phi }_1\omega _0^2R_{}\varsigma \right)\left[\frac{1+2\kappa _{\mathrm{eff}}}{2(1+\kappa _{\mathrm{eff}})}\right]+\frac{m}{G_{22}n_2}\left(i\omega \mathrm{\Phi }_2\omega _0^2R_{}\varsigma \right)\left[\frac{2+\kappa _{\mathrm{eff}}}{2(1+\kappa _{\mathrm{eff}})}\right]\right),$$

where $`\sigma _0`$ is the tension for the stationary configuration. As a result, the boundary condition (18) contains extra terms and becomes

$$\left(mn_1\frac{2m\sigma _0}{R_{}G_{11}n_1}\left[\frac{1+2\kappa _{\mathrm{eff}}}{2(1+\kappa _{\mathrm{eff}})}\right]\right)\left[\omega _0^2R_{}\frac{\mathrm{\Phi }_1}{r}\omega ^2\mathrm{\Phi }_1\right]\frac{\sigma _0}{R_{}^2}\left(2\frac{\mathrm{\Phi }_1}{r}+\mathrm{\Delta }_{\theta ,\varphi }\frac{\mathrm{\Phi }_1}{r}\right)=$$

$$=\left(mn_2+\frac{2m\sigma _0}{R_{}G_{22}n_2}\left[\frac{2+\kappa _{\mathrm{eff}}}{2(1+\kappa _{\mathrm{eff}})}\right]\right)\left[\omega _0^2R_{}\frac{\mathrm{\Phi }_2}{r}\omega ^2\mathrm{\Phi }_2\right],$$
(27)
while Eq. (19) remains the same. For the radial part of the velocity potential we obtain

$$\left(1\frac{\sigma _0}{R_{}P}\left[\frac{1+2\kappa _{\mathrm{eff}}}{2(1+\kappa _{\mathrm{eff}})}\right]\right)\left[\omega _0^2R_{}\frac{\mathrm{\Phi }_1}{r}\omega ^2\mathrm{\Phi }_1\right]+\frac{\sigma _0(l1)(l+2)}{mn_1R_{}^2}\frac{\mathrm{\Phi }_1}{r}=$$

$$=\kappa _{\mathrm{eff}}\left(1+\frac{\sigma _0}{R_{}P}\left[\frac{2+\kappa _{\mathrm{eff}}}{2(1+\kappa _{\mathrm{eff}})}\right]\right)\left[\omega _0^2R_{}\frac{\mathrm{\Phi }_2}{r}\omega ^2\mathrm{\Phi }_2\right],$$
(28)
where $`P`$ is the pressure at the interface. The modified boundary condition results in the following equation for the low frequency modes

$$\omega ^2\omega _0^2\left(1\frac{3\sigma _0}{2R_{}P}+\frac{\sigma _0(l1)(l+2)}{mn_1\omega _0^2R_{}^3}\kappa _{\mathrm{eff}}\right)f.$$
(29)
where $`n_1`$ is the density at the interface. The interface tension shifts the frequencies of the lowest modes and narrows the stability region. For $`l=1`$ the inner droplet moves as a whole without changing its shape. However, the displacement of the droplet into the less dense region decreases the interface energy which is proportional to $`n^{3/2}`$. This is the origin for the contribution $`3\sigma _0/R_{}P`$ to the mode frequency. The mode frequencies now become imaginary at

$$\kappa ^2>1\frac{\sigma _0}{R_{}P}+\frac{2\sigma _0(l1)(l+2)}{mn_1\omega _0^2R_{}^3}.$$
(30)
As $`\kappa `$ increases the mode with $`l=1`$ becomes imaginary first which determines the system’s stability limit. When the surface tension is absent, the higher $`l`$ is the more unstable the mode is. In the present case, the two component condensate is locally unstable when

$$\frac{a_{11}}{a_{22}}>1\frac{\sigma _0}{R_{}P}.$$
(31)
## V Global stability: Energy of asymmetric phase

So far we have discussed the local stability of the symmetric phase. We next turn our attention to its global stability by comparing the energies of the symmetric phase and the asymmetric phase. Let us consider an asymmetric state in which component 1 sits on the top of the component 2. The transition occurs when $`\kappa `$ is close to $`1`$ which we assume. Hence, the interface is approximately flat (see Fig. 8).

Fig. 8. Symmetric and Asymmetric phases of two component BECs.

In the Thomas-Fermi approximation, the density distribution is still given by Eqs. (4) and (5), but the interface position is now different. The TF density distribution suggests that the outer edge of the condensates has a spherical shape and, therefore, in the asymmetric configuration the condensates have equal chemical potentials $`\mu _{1a}=\mu _{2a}=\mu `$. Carrying out the integration, we find that the energy of the asymmetric state is given by

$$E_a=\frac{\mu }{2}\left(N_1+N_2\right)+\frac{\pi m^2\omega _0^4R^7}{8}\left(\frac{1}{G_{11}}\frac{1}{G_{22}}\right)\left[\frac{4}{35}\frac{d}{6R}+\frac{d^5}{10R^5}\frac{d^7}{21R^7}\right]+\frac{\pi m^2\omega _0^4R^7}{35G_{22}}+E_{sa},$$
(32)
where $`z=d`$ is the position of the interface, $`E_{sa}`$ is the interface energy and $`R=\sqrt{2\mu /m\omega _0^2}`$ is the radius of the sphere. The radius $`R`$ (and hence $`\mu `$) and the position $`d`$ are determined by the number of particles

$$N_1=\frac{\pi \mu R^3}{15G_{11}}\left[4\frac{15d}{2R}+\frac{5d^3}{R^3}\frac{3d^5}{2R^5}\right],$$
(33)
$$N_2=\frac{\pi \mu R^3}{15G_{22}}\left[4+\frac{15d}{2R}\frac{5d^3}{R^3}+\frac{3d^5}{2R^5}\right].$$
(34)
From Eqs. (33), (34), we obtain

$$G_{11}N_1+G_{22}N_2=\frac{4\pi m\omega _0^2R^5}{15},$$
(35)
which determines $`R`$. The position of the interface can be found from

$$4\frac{15d}{2R}+\frac{5d^3}{R^3}\frac{3d^5}{2R^5}=\frac{8G_{11}N_1}{G_{22}N_2+G_{11}N_1}.$$
(36)
The energy of a spherically symmetric configuration when the first droplet is located inside and the second outside is given by

$$E_0=\frac{\mu _1N_1}{2}+\frac{\mu _2N_2}{2}+\pi m\omega _0^2[(\frac{\mu _1}{G_{11}}\frac{\mu _2}{G_{22}})\frac{R_{}^5}{5}+\frac{m\omega _0^2}{2}(\frac{1}{G_{22}}\frac{1}{G_{11}})\frac{R_{}^7}{7}+$$

$$+\frac{2^{7/2}\mu _2^{7/2}}{35G_{22}m^{5/2}\omega _0^5}]+E_{s0},$$
(37)
where $`E_{s0}`$ is the interface energy in the symmetric state.

When the system is close to the point of phase transition the interface energies of two configurations are given by (for details see Sec. VI)

$$E_{s0}=\frac{8\pi \mathrm{}R_{}^2\mu ^{3/2}}{\sqrt{3m}G}\sqrt{a_{12}/\sqrt{a_{11}a_{22}}1}\left(1\frac{R_{}^2}{R^2}\right)^{3/2},$$
(38)
$$E_{sa}=\frac{4\pi \mathrm{}R^2\mu ^{3/2}}{5\sqrt{3m}G}\sqrt{a_{12}/\sqrt{a_{11}a_{22}}1}\left(1\frac{d^2}{R^2}\right)^{5/2}.$$
(39)
In general, these expressions need to be numerically evaluated. They are analytically tractable in the limit when one of the components is particularly small, to which we turn our attention next.

### A Energy for $`N_1N_2`$

If the number of atoms of the first species is much less than $`N_2`$ then $`R_{}R`$, $`d/R1(4N_1/5N_2)^{1/3}`$ and from Eqs. (38), (39) we find

$$E_{s0}\frac{8\pi \mathrm{}R_{}^2\mu ^{3/2}}{\sqrt{3m}G}\sqrt{a_{12}/\sqrt{a_{11}a_{22}}1},$$
(40)
$$E_{sa}1.008\sqrt{\frac{R_{}}{R}}E_{s0}E_{s0}.$$

The interface energy of the asymmetric state is negligible. Close to the point of phase transition the energy difference between the asymmetric and symmetric configuration is given by

$$E_aE_0=\frac{1}{2}\mu N_1\left[(1\kappa )^2\frac{4\sqrt{6}\xi }{R_{}}\sqrt{a_{12}/\sqrt{a_{11}a_{22}}1}\right]\text{,}$$
(41)
where $`\xi =\mathrm{}/\sqrt{2m\mu }`$.

In the absence of the interface tension ($`\xi 0`$) the asymmetric configuration has greater energy no matter what is the value of $`\kappa `$ in comparison with $`1`$. However, the interface tension makes the asymmetric phase preferable at some $`\kappa `$ and results in symmetric-asymmetric phase transition. The condition $`E_a=E_0`$ describes a line of global instability of the system

$$\frac{a_{11}}{a_{22}}=14\left(\frac{\sqrt{6}\xi }{R_{}}\right)^{1/2}(a_{12}/\sqrt{a_{11}a_{22}}1)^{1/4}\text{.}$$
(42)
The soft mode condition gives the line of local stability

$$\frac{a_{11}}{a_{22}}=1\frac{4\sqrt{6}\xi }{3R_{}}(a_{12}/\sqrt{a_{11}a_{22}}1)^{1/2}\text{.}$$
(43)
In Fig. 9 we plot the phase diagram that shows different stability regions of two component condensates. In estimates we take $`\xi /R_{}=0.01`$. When $`a_{12}<\sqrt{a_{11}a_{22}}`$ the homogeneous binary mixture is a stable state. Otherwise the two components are phase separated. In the later case for small ratio $`a_{11}/a_{22}`$ the state $`(1,2)`$ with the first component being inside and the second outside is the only stable configuration. If the ratio $`a_{11}/a_{22}`$ increases the configuration becomes globally unstable when we cross the left (solid) curve. However, the system is locally stable since the configuration corresponds to a local minima of energy. Further increase of $`a_{11}/a_{22}`$ crosses the line of local stability (another solid line) and the system undergoes a phase transition into a new stable asymmetric state. If initially the system is prepared in the $`(2,1)`$ state then with decreasing the ratio $`a_{11}/a_{22}`$ the configuration first becomes globally unstable when we cross the right dotted line and locally unstable when $`a_{11}/a_{22}`$ is close enough to $`1`$.

Fig. 9. Phase diagram of the two component condensate in coordinates displaying relative interaction strength between bosons.

## VI Energy as a function of droplet displacement

The question about the system’s global stability is related to the energy of the state. To further explore the metastability of the system here we estimate the energy of the system for a displaced inner droplet. We assume the inner droplet has a spherical shape with a radius $`r_0`$ and the center of the sphere is displaced a distance $`z_0`$ from the trap center. The radius $`r_0`$ depends on $`z_0`$. The condensate densities are given by the same TF expressions (4), (5), however, the position of the interface is now different. We assume the number of particles of each species are fixed:

$$N_1=\frac{4\pi }{3G_{11}}\left[\mu _1r_0^3\frac{3m\omega _0^2r_0^5}{10}\frac{m\omega _0^2r_0^3z_0^2}{2}\right],$$
(44)
$$N_2=\frac{4\pi }{3}\frac{1}{G_{22}}\left[\frac{2}{5}\mu _2R_2^3\mu _2r_0^3+\frac{3m\omega _0^2r_0^5}{10}+\frac{m\omega _0^2r_0^3z_0^2}{2}\right].$$
(45)
For fixed numbers of particles, the chemical potentials are functions of the positions of the droplets and are not constants. In the TF limit we find for the energy of the system

$$E=\frac{\mu _1N_1}{2}+\frac{\mu _2N_2}{2}+\pi m\omega _0^2[(\frac{\mu _1}{G_{11}}\frac{\mu _2}{G_{22}})(\frac{r_0^5}{5}+\frac{r_0^3z_0^2}{3})+\frac{m\omega _0^2}{2}(\frac{1}{G_{22}}\frac{1}{G_{11}})[\frac{r_0^7}{7}+\frac{2r_0^5z_0^2}{3}+\frac{r_0^3z_0^4}{3}]+$$

$$+\frac{2^{7/2}\mu _2^{7/2}}{35G_{22}m^{5/2}\omega _0^5}]+E_s,$$
(46)
where $`E_s`$ is the interface energy. The energy (46) depends on four parameters: $`\mu _1`$, $`\mu _2`$, $`z_0`$, $`r_0`$. However, only two parameters are independent because there are two restrictions imposed by particle conservation (44), (45). For simplicity let us consider the path for which $`\mu _2=const`$. Then only one parameter is independent. If we choose $`r_0`$ as an independent variable then the energy as a function of $`r_0`$ can be expressed in an explicit form

$$E=const+\frac{2\pi (\mu _{10}\mu _2)^2R_{}^6}{3G_{11}r_0^3}+$$

$$+4\pi \mu _2^2\left(\frac{1}{G_{22}}\frac{1}{G_{11}}\right)\left[\frac{R_{}^6}{6r_0^3}+\frac{2r_0^2(r_0^3R_{}^3)}{15R_2^2}\frac{R_{}^8}{5R_2^2r_0^3}+\frac{3R_{}^{10}}{50R_2^4r_0^3}+\frac{2r_0^2R_{}^5}{25R_2^4}\frac{12r_0^7}{175R_2^4}\right]+E_s,$$
(47)
where $`\mu _2=m\omega _0^2R_2^2/2`$, $`\mu _{10}`$ is the chemical potential at $`z_0=0`$ and $`R_{}`$ is the droplet radius at $`z_0=0`$. The dependence of the droplet radius on displacement $`z_0`$ is described by Eq. (45).

### A Interface energy

Here we estimate the interface energy as a function of the droplet displacement $`z_0`$. According to Ref.  the interface tension is given by

$$\sigma =\frac{4}{\sqrt{3}}\sqrt{(\xi _1^2+\xi _2^2)\left[a_{12}/\sqrt{a_{11}a_{22}}1\right]}P,$$
(48)
where $`\xi _i`$ represents the single condensate coherence length $`\xi _i=\mathrm{}/\sqrt{2m_iG_{ii}n_i}`$, the pressure $`PG_{ii}n_i^2/2`$ and the condensate densities $`n_i`$ are estimated near the interface. The interface tension $`\sigma `$ depends on the local density and, hence, changes along the interface. To calculate $`\sigma `$ we assume that the system is close to the point of phase transition, that is $`G_{11}G_{22\text{ }}=G`$, $`R_1R_2=R`$, $`\mu _1\mu _2=\mu `$ and near the interface $`n_1n_2=n`$. As a result,

$$\sigma \frac{2\mathrm{}\sqrt{G}n^{3/2}}{\sqrt{3}\sqrt{m}}\sqrt{a_{12}/\sqrt{a_{11}a_{22}}1}$$
(49)
and we find for the surface energy

$$E_s=\frac{4\pi \mathrm{}\mu ^{3/2}R^2r_0}{5\sqrt{3m}Gz_0}\sqrt{a_{12}/\sqrt{a_{11}a_{22}}1}\left[\left(1\frac{(r_0z_0)^2}{R^2}\right)^{5/2}\left(1\frac{(r_0+z_0)^2}{R^2}\right)^{5/2}\right].$$
(50)
### B Energy for small inner droplet

In the limit $`r_0R`$ ($`N_1N_2`$) Eq. (45) results in a simple relation between the droplet radius $`r_0`$ and its displacement $`z_0`$: $`r_0^3\mu _2R_{}^3/(\mu _2m\omega _0^2z_0^2/2)`$. Then near the point of phase transition the energy of the system is given by

$$Econst+N_1\mu \left[(1\kappa )\frac{z_0^2}{R^2}+b\left(1\frac{z_0^2}{R^2}\right)^{5/6}\right],$$
(51)
where

$$b=\frac{2\sqrt{6}\xi }{R_{}}\sqrt{a_{12}/\sqrt{a_{11}a_{22}}1}=\frac{3\sigma _0}{2PR_{}}.$$

The first term in Eq. (51) corresponds to energy contribution due to difference between $`a_{11}`$ and $`a_{22}`$, it increases with $`z_0`$ when $`\kappa <1`$. The second term is the contribution from the interface energy, it decreases with increasing $`z_0`$ because surface tension drops when the droplet moves in the region with smaller density. In Fig. 10 we plot the energy (51) as a function of displacement of the inner droplet $`z_0`$. The energy is shown for $`b=0.3`$ and different $`\kappa `$. Far from the point of phase transition the energy monotonically increases with displacement of the droplet, in this case the first term dominates. When $`\kappa `$ goes closer to $`1`$ the first term in Eq. (51) becomes small and the interface contribution changes the energy behavior. The energy becomes nonmonotonic function of displacement. The value of $`\kappa `$ when the final and initial energies are equal is the onset of global instability. However, the system is still locally stable since position of the droplet at the center corresponds to a local minima of energy. In this region an energy barrier prevents the droplet to move outside. Further increase of $`\kappa `$ results in disappearance of the energy barrier and initial configuration becomes both locally and globally unstable.

Fig. 10. Energy of the system as a function of displacement of the inner droplet for fixed $`b=0.3`$ and different $`\kappa `$.

Near the point of local instability the energy barrier is $`\mathrm{\Delta }E18N_1\mu (\kappa _s\kappa )^2/5b`$, where $`\kappa _s=15b/6`$. The barrier becomes of the order of $`\mu `$ when $`\kappa _s\kappa \sqrt{5b/18N_1}`$. If $`N_110^6`$, $`b0.1`$ we obtain $`\kappa _s\kappa 10^4`$, while $`1\kappa _s0.1`$.

## VII Normal modes for different trap frequencies and particle masses of the two species

In this section we generalize our results to the case $`\omega _1\omega _2`$, $`m_1m_2`$ which is useful for future experiments where the trap frequencies and particles can differ for each component. For simplicity we omit the interface tension and consider spherical trapping potentials. The generalized boundary condition for the velocity potential $`\mathrm{\Phi }`$ reads

$$m_1\omega _1^2R_{}\frac{\mathrm{\Phi }_1}{r}m_1\omega ^2\mathrm{\Phi }_1=\kappa m_2\omega _2^2R_{}\frac{\mathrm{\Phi }_2}{r}\kappa m_2\omega ^2\mathrm{\Phi }_2,$$
(52)
where $`\kappa =\sqrt{G_{11}/G_{22}}=\sqrt{a_{11}m_2/a_{22}m_1}`$. The boundary condition for $`\mathrm{\Phi }_1/r`$ remains the same as Eq. (19). Now the phase transition is controlled by a new dimensionless parameter

$$\stackrel{~}{\kappa }=\frac{m_2\omega _2^2}{m_1\omega _1^2}\sqrt{\frac{G_{11}}{G_{22}}}=\frac{m_2^{3/2}\omega _2^2}{m_1^{3/2}\omega _1^2}\sqrt{\frac{a_{11}}{a_{22}}}.$$

Instead of (23), the solutions for the inner and outer condensates have the form

$$\mathrm{\Phi }_1=C_1r^lF(\alpha _1,\beta _1,l+3/2,\frac{r^2}{R_1^2}),\mathrm{\Phi }_2=C_2r^lF(\alpha _2,\beta _2,1,1\frac{r^2}{R_2^2}),$$
(53)
where $`R_i^2=2\mu _i/m_i\omega _i^2`$ ($`i=1,2`$) and

$$\alpha _i=\frac{1}{2}\left[l+\frac{3}{2}\sqrt{l^2+l+\frac{9}{4}+\frac{2\omega ^2}{\omega _i^2}}\right],\beta _i=\frac{1}{2}\left[l+\frac{3}{2}+\sqrt{l^2+l+\frac{9}{4}+\frac{2\omega ^2}{\omega _i^2}}\right].$$

Using Eq. (53) and the boundary conditions we obtain an equation for the normal mode frequencies

$$\omega ^2=\frac{\omega _1^2(1\stackrel{~}{\kappa })\left[l(l+3/2)+(l\omega ^2/\omega _1^2)\lambda xs_1(\omega ,x)\right]\left[(l\omega ^2/\omega _2^2)xs_2(\omega ,x)l\right]}{\left[\stackrel{~}{\kappa }(l\omega ^2/\omega _1^2)\lambda xs_1(\omega ,x)\omega _1^2/\omega _2^2+(l+3/2)\left(l(\stackrel{~}{\kappa }\omega _1^2/\omega _2^21)+(l\omega ^2/\omega _2^2)xs_2(\omega ,x)\right)\right]},$$
(54)
where

$$\lambda =\frac{m_1\omega _1^2\mu _2}{m_2\omega _2^2\mu _1},x=\frac{R_{}^2}{R_2^2},s_1(\omega ,x)=\frac{F(\alpha _1+1,\beta _1+1,l+5/2,\lambda x)}{F(\alpha _1,\beta _1,l+3/2,\lambda x)},s_2(\omega ,x)=\frac{F(\alpha _2+1,\beta _2+1,2,1x)}{F(\alpha _2,\beta _2,1,1x)}.$$

In terms of $`\stackrel{~}{\kappa }`$ and $`\lambda `$ the position of the interface is given by

$$R_{}=R_2\sqrt{\frac{1\stackrel{~}{\kappa }\lambda }{\lambda (1\stackrel{~}{\kappa })}}.$$

In the region $`|1\stackrel{~}{\kappa }|1`$one can take $`\omega =0`$ in the right side of Eq. (54) and put $`\stackrel{~}{\kappa }`$,$`\lambda 1`$ in the multiple. As a result, we find

$$\omega ^2\omega _1^2(1\stackrel{~}{\kappa })f(l,x),$$
(55)
where

$$f(l,x)=\frac{l\left(l+3/2+xs_1(0,x)\right)\left(xs_2(0,x)1\right)}{\left[xs_1(0,x)\omega _1^2/\omega _2^2+(l+3/2)\left(\omega _1^2/\omega _2^21+xs_2(0,x)\right)\right]}>0.$$

Eq. (55) is a generalization of Eq. (25) for the low frequency modes. The point of instability is now determined by the condition $`\stackrel{~}{\kappa }=1`$, or $`a_{11}m_2^3\omega _2^4=a_{22}m_1^3\omega _1^4`$. $`\stackrel{~}{\kappa }`$ can be changed not only by changing the relative interaction strengths, but also by changing the ratio of the trapping frequencies. Furthermore it is proportional to a higher (4th) power of the ratio $`\omega _1/\omega _2`$ as compared to $`a_{11}/a_{22}`$ and thus will be very sensitive to the frequency change. Near the instability point $`\omega \sqrt{1\stackrel{~}{\kappa }}`$ and becomes imaginary when $`\stackrel{~}{\kappa }>1`$.

In summary, we study the normal modes and the stability of the phase separated two component BECs and predict a new symmetric-asymmetric transition in such a system. We show that the normal mode frequencies differ from those for the one component condensate and new branches appear. We found normal modes with low frequencies which are analogous to gravitational waves at the interface between two immiscible fluids. Under the change of the relative interaction strength $`a_{11}/a_{22}`$ or the ratio $`\omega _1/\omega _2`$ the frequencies of those modes go to zero and become imaginary. Imaginary frequencies mean instability and the system undergoes a quantum phase transition into a new stable configuration. Such transition can be observed experimentally by manipulating the interaction strength via Feschbach resonances or by changing the trap frequencies. Interface tension shifts the normal mode frequencies and results in stability of the asymmetric state in some parameter region. Another interesting possibility is to study tunneling by turning on a microwave field which allows direct exchange of the two components and has been done in many JILA experiments with two hyperfine states of 87Rb. If we suddenly switch the species, one can study a decay of new unstable configuration or tunneling directly, instead of trying to tune the ratio $`a_{11}/a_{22}`$.

We are grateful to A. Fetter for valuable remarks. This work was supported by NASA, Grant No. NAG8-1427.


text:
# Theoretical description of the ferromagnetic 𝜋-junctions near the critical temperature

## I Introduction

The strong exchange field acting on the electrons in a ferromagnet provokes the damped oscillatory behavior of the superconducting order parameter. This effect is at the origin of the $`\pi `$-junction realization in superconductor-ferromagnet-superconductor (S/F/S) systems . The study of the superconducting $`\pi `$-state sheds more light on the coexistence of superconductivity and magnetism in general and may be also important for superconducting electronics (see for example ).

Recently such $`\pi `$-junctions have been successfully fabricated and experimentally studied in . While the existing theoretical approach  provides a good qualitative description of the observed damped oscillations of the critical current as a function of ferromagnetic layer thickness , a more detailed description is needed to take into account the finite S/F boundary transparency. Indeed, in  Nb/Al/Al2O3/PdNi/Nb junctions have been studied. In such junctions, the presence of the Al2O3 barrier can be modeled by a very low transparency of one of the S/F interfaces.

In this work we develop a theory of ferromagnetic junctions with different transparencies of the interfaces providing a general description of such junctions near the critical temperature.

## II General formalism

We will concentrate on the studies of the properties of an S/F/S junction with a thin F-layer of thickness $`d`$ and large superconducting electrodes, see Fig. $`1`$. The very convenient set of equations describing an inhomogeneous superconductivity has been elaborated by Eilenberger . These are transport-like equations for the energy-integrated Green’s functions $`f`$ and $`g`$, assuming that relevant length scales are much larger than atomic length scales. For our geometry all quantities depend only on one coordinate $`x`$, and in the presence of the ferromagnetic exchange field $`h(x)`$ acting on the electron spins in the F-layer, the Eilenberger equations take the form (see for example )

$$\left(\omega +ih(x)+\frac{1}{2\tau }G(x,\omega )\right)f(x,\theta ,\omega )+\frac{1}{2}v_F\mathrm{cos}\theta \frac{f(x,\theta ,\omega )}{x}=\left(\mathrm{\Delta }(x)+\frac{1}{2\tau }F(x,\omega )\right)g(x,\theta ,\omega ),$$
(1)
$$G(x,\omega )=\frac{d\mathrm{\Omega }}{4\pi }g(x,\theta ,\omega ),F(x,\omega )=\frac{d\mathrm{\Omega }}{4\pi }f(x,\theta ,\omega ),$$
(2)
$$f(x,\theta ,\omega )f^+(x,\theta ,\omega )+g^2(x,\theta ,\omega )=1,$$
(3)
where $`\omega =2\pi T\left(n+1/2\right)`$ are the Matsubara frequencies and $`\theta `$ is the angle between the $`𝐱`$ axis and the direction of the Fermi velocity $`v_F`$, and $`\tau `$ is the elastic scattering time. In the following, we consider the behavior of S/F/S systems close to the critical transition temperature $`T_c`$, so we may linearize equations $`\left(\text{2}\right)`$ by putting $`g=sign\left(\omega \right).`$

Usually the electron scattering mean free path in S/F/S systems is rather small. In such a dirty limit the angular dependence of the Green’s functions is weak, and the Eilenberger equations can be replaced by much simpler Usadel equations . In fact, the conditions of the applicability of Usadel equations are $`T_c\tau 1`$ and $`h\tau 1.`$ The second condition, is much more restrictive due to a large value of the exchange field ($`hT_c`$). Therefore the attempts to retain in the Usadel equations the first correction in the parameter $`h\tau `$ have been made in . They have resulted in the renormalization of the diffusion constant in F-layer $`D_fD_f(12ih\tau ),`$ where $`D_f=v_F^2\tau /3.`$ The critical analysis of this renormalization in  has revealed its inaccuracy, but has not provided the right correction. To resolve this controversy we perform below a careful derivation of the Usadel equation for an F-layer retaining the linear correction over parameter $`h\tau .`$

In the ferromagnet the Cooper pairing is absent, so $`\mathrm{\Delta }=0`$, and considering the case of the strong exchange field$`\omega h`$, we may write the linearized Eilenberger equation for a function $`f`$ as

$$\left(ih+\frac{1}{2\tau }sign\left(\omega \right)\right)f+\frac{1}{2}v_F\mathrm{cos}\theta \frac{f}{x}=\frac{1}{2\tau }sign\left(\omega \right)F,$$
(4)
where $`v_F`$ and $`\tau `$ refer to the Fermi velocity and the scattering time in the ferromagnet. Let us recall that $`F`$ is an anomalous Green function averaged over the angles. Then we may present $`f`$ as $`f=F(x)+f_1(x,\theta ),`$ where $`f_1F`$. So for $`\omega >0`$

$$ihf+\frac{f_1}{2\tau }+\frac{1}{2}v_F\mathrm{cos}\theta \frac{f}{x}=0.$$
(5)
The averaging of this equation over the angles gives

$$ihF+\frac{1}{2}v_F\mathrm{cos}\theta \frac{f_1}{x}=0.$$
(6)
Multiplying $`\left(\text{5}\right)`$ by cos$`\theta `$, averaging and taking the derivative over $`x`$ we obtain

$$\mathrm{cos}\theta \frac{f_1}{x}=\frac{\tau v_F}{\left(2ih\tau +1\right)}\left(\mathrm{cos}^2\theta \frac{^2f_1}{x^2}+\frac{1}{3}\frac{^2F}{x^2}\right).$$
(7)
Then combining $`\left(\text{5}\right)`$ and $`\left(\text{6}\right)`$ we have

$$ihF\frac{\tau v^2}{2\left(2ih\tau +1\right)}\left(\mathrm{cos}^2\theta \frac{^2f_1}{x^2}+\frac{1}{3}\frac{^2F}{x^2}\right)=0.$$
(8)
Neglecting in this equation the $`f_1`$ term, we obtain the usual Usadel equation with the renormalized diffusion constant $`D_fD_f(12ih\tau ).`$ However the omitted term gives a correction of the same order as that coming from the diffusion constant renormalization . It effectively means that it is needed to keep higher order terms in the expansion of $`f`$ over $`\mathrm{cos}\theta `$ (Legendre polynomials). Multiplying $`\left(\text{5}\right)`$ consecutively by $`\mathrm{cos}^2\theta `$ and $`\mathrm{cos}^3\theta `$ and taking derivative over $`x`$, we obtain after some algebra the following exact equation

$$ihF\frac{\tau v_F^2}{2\left(2ih\tau +1\right)}\left\{\frac{1}{\left(2ih\tau +1\right)}\left[\frac{2ih\tau }{3}\frac{^2F}{x^2}\frac{v_F^2\tau ^2}{\left(2ih\tau +1\right)}\left(\mathrm{cos}^4\theta \frac{^4f_1}{x^4}+\frac{1}{5}\frac{^4F}{x^4}\right)\right]+\frac{1}{3}\frac{^2F}{x^2}\right\}=0.$$
(9)
Here we may already safely neglect the $`f_1`$ term, as it gives a contribution $`\left(h\tau \right)^2`$. In the result, in the linear approximation over $`h\tau ,`$ the Usadel equation is written as

$$ihF\frac{\tau v_F^2}{6\left(4ih\tau +1\right)}\frac{^2F}{x^2}\frac{\tau ^3v_F^4}{10}\frac{^4F}{x^4}=0.$$
(10)
Taking in mind that in zero approximation over $`h\tau ,`$ the function $`F`$ satisfies the equation

$$ihF\frac{\tau v_F^2}{6}\frac{^2F}{x^2}=0,$$
(11)
we may finally rewrite $`\left(\text{10}\right)`$ in the following form applicable for positive and negative $`\omega `$

$$isign\left(\omega \right)hF\frac{D_f\left[1i(2/5)h\tau sign\left(\omega \right)\right]}{2}\frac{^2F}{x^2}=0.$$
(12)
So we conclude that the first correction in the parameter $`h\tau `$ to the Usadel equation leads to the somewhat different renormalization of the diffusion constant $`D_fD_f(1i(2/5)h\tau sign\left(\omega \right))`$, comparing to what has been suggested before  $`\left(D_fD_f\left[1i2h\tau sign\left(\omega \right)\right]\right)`$. This result has been also independently obtained by Tagirov .

In the superconducting layer, the linearized Usadel equation is

$$\left|\omega \right|F(x,\omega )\frac{D_s}{2}\frac{^2F(x,\omega )}{x^2}=\mathrm{\Delta }(x),$$
(13)
$$\mathrm{\Delta }(x)=\left|\lambda \right|\pi T\underset{\omega }{}F(x,\omega ),$$
(14)
where $`\lambda `$ is the BCS coupling constant and $`D_s`$ is the diffusion coefficient in S layer.

The equations $`\left(\text{12}\right)`$ and $`\left(\text{13}\right)`$ are completed by the general boundaries conditions at the S/F interfaces

$`F_s\left(d/2\right)=F_f\left(d/2\right)+\xi _f\gamma _{B2}\left({\displaystyle \frac{F_f}{x}}\right)_{d/2},`$

$$F_s\left(d/2\right)=F_f\left(d/2\right)\xi _f\gamma _{B1}\left(\frac{F_f}{x}\right)_{d/2},$$
(15)
$`\left({\displaystyle \frac{F_s}{x}}\right)_{\pm d/2}={\displaystyle \frac{\sigma _n}{\sigma _s}}\left({\displaystyle \frac{F_f}{x}}\right)_{\pm d/2},`$

where the notation $`F_s(F_f)`$ is used for the anomalous Green function in a superconductor (ferromagnet) and $`\sigma _n`$ $`\left(\sigma _s\right)`$ is the conductivity of the F-layer $`\left(\text{S-layer above }T_c\right),`$ $`\xi _f=\sqrt{\frac{D_f}{h}\text{ }}`$is the characteristic length of superconducting correlations decaying in F-layer, and $`\xi _s=\sqrt{\frac{D_s}{2T_c}}`$ is the superconducting coherence length of the S-layer, the parameter $`\gamma _B=\frac{R_b\sigma _f}{\xi _f},`$ where $`R_b`$ is the S/F boundary resistance per unit area. Here we introduce two different parameter $`\gamma _{B1}`$ and $`\gamma _{B2}`$, for the left and right S/F interfaces, assuming that their boundary resistances may be different. Note that the parameter $`\gamma _B`$ is directly related to the transparency of the S/F interface $`T=\frac{1}{1+\gamma _B}`$ . The limit $`T=0`$ $`\left(\gamma _B=\mathrm{}\right)`$ corresponds to a vanishingly small boundary transparency, and the limit $`T=1`$ $`\left(\gamma _B=0\right)`$ corresponds to a perfectly transparent interface.

In principle, the equations $`\left(\text{12}\right)`$ and $`\left(\text{13}\right)`$ with boundary conditions $`\left(\text{15}\right)`$ give the complete description of the S/F/S junction near the transition temperature.

## III Calculation of the Josephson critical current

Let us start with a general solution of the Usadel equation in the ferromagnet for $`\omega >0`$

$$F_f\left(\omega >0,x\right)=A_\omega \mathrm{sinh}(kx)+B_\omega \mathrm{cosh}(kx),$$
(16)
where $`k=[(10.2h\tau )+i(1+0.2h\tau ]\sqrt{\frac{h}{D_f}}.`$ Similarly for $`\omega <0`$

$$F_f\left(\omega <0,x\right)=C_\omega \mathrm{sinh}(k^{}x)+D_\omega \mathrm{cosh}(k^{}x)$$
(17)
We see that the complex wave-vector $`k`$ describes the oscillating exponential damping of the anomalous Green function inside the ferromagnet. In the case of very small scattering time $`\left(\tau 0\right),`$ the imaginary and real parts of the wave-vector $`k`$ coincide, and so the characteristic lengths for oscillating and damping are exactly the same. The finite scattering length decreases the oscillating period and increases the damping length. This is consistent with the fact that in a pure limit the exponential damping is replaced by a much slower power damping, while the oscillations have a much shorter period comparing to the dirty limit .

To obtain the analytical description of the critical current in S/F/S junctions, it is useful to introduce the functions $`F^+`$ and $`F^{}`$

$`F^+\left(\omega >0,x\right)`$ $`=`$ $`F\left(\omega >0,x\right)+F\left(\omega <0,x\right),`$ (18)
$`F^{}\left(\omega >0,x\right)`$ $`=`$ $`F\left(\omega >0,x\right)F\left(\omega <0,x\right).`$ (19)
In a superconductor, the equation for $`F_s^+`$ is the same as $`\left(\text{13}\right)`$ with $`\mathrm{\Delta }2\mathrm{\Delta }`$, while the equation for $`F_s^{}`$ is much simpler

$$\left|\omega \right|F_s^{}(x,\omega )\frac{D_s}{2}\frac{^2F_s^{}(x,\omega )}{x^2}=0,$$
(20)
and the self consistency equation reads as

$$\mathrm{\Delta }=\left|\lambda \right|\pi T\underset{\omega >0}{}F_s^+(x,\omega ).$$
(21)
The boundary conditions $`\left(\text{15}\right)`$ will be the same for $`F^+`$ and $`F^{}`$ functions too. Using the solution of $`\left(\text{20}\right)`$ for $`F_s^{}(x,\omega )`$, it may be demonstrated (see also) that in the case of a rather large resistivity of the F layer $`\left(\frac{\sigma _n}{\sigma _s}\frac{\xi s}{\xi _f}1\right)`$ or in the case low transparency of the interfaces, when the condition $`\frac{\sigma _n}{\sigma _s}\frac{\xi s}{\xi _f}\frac{1}{\gamma __B}1`$ is satisfied, the boundary conditions for $`F_f^{}`$ functions are essentially simplified

$`F_f^{}\left(d/2\right)+\xi _f\gamma _{B2}\left({\displaystyle \frac{F_f^{}}{x}}\right)_{_{d/2}}`$ $`=`$ $`0,`$ (22)
$`F_f^{}\left(d/2\right)\xi _f\gamma _{B1}\left({\displaystyle \frac{F_f^{}}{x}}\right)_{_{d/2}}`$ $`=`$ $`0.`$ (23)
Further on, introducing the functions

$`^+`$ $`=`$ $`\left|\lambda \right|\pi T{\displaystyle \underset{\omega >0}{}}F^+\left(\omega >0,x\right),`$ (24)
$`^{}`$ $`=`$ $`\left|\lambda \right|\pi T{\displaystyle \underset{\omega >0}{}}F^{}\left(\omega >0,x\right),`$ (25)
we obtain the following system of equations for the order parameter and its derivatives on the left and right sides of the junction

$$_f^{}\left(d/2\right)+\xi _f\gamma _{B2}\left(\frac{_f^{}}{x}\right)_{d/2}=0,$$
(26)
$$_f^{}\left(d/2\right)\xi _f\gamma _{B1}\left(\frac{_f^{}}{x}\right)_{d/2}=0,$$
(27)
$$\mathrm{\Delta }\left(x=d/2\right)=\mathrm{\Delta }^+=_f^+\left(d/2\right)+\xi _f\gamma _{B2}\left(\frac{_f^+}{x}\right)_{d/2},$$
(28)
$$\mathrm{\Delta }\left(x=d/2\right)=\mathrm{\Delta }^{}=_f^+\left(d/2\right)\xi _f\gamma _{B1}\left(\frac{_f^+}{x}\right)_{d/2},$$
(29)
$$\left(\frac{\mathrm{\Delta }^+}{x}\right)=\frac{\sigma _n}{\sigma _s}\left(\frac{_f^+}{x}\right)_{d/2},$$
(30)
$$\left(\frac{\mathrm{\Delta }^{}}{x}\right)=\frac{\sigma _n}{\sigma _s}\left(\frac{_f^+}{x}\right)_{d/2}.$$
(31)
This system of six equations, after the elimination of four coefficients of the type $`\left(\left|\lambda \right|\pi T\underset{\omega >0}{}A_\omega \right)`$ may be reduced to a standard linear form of $`\left(\frac{\mathrm{\Delta }^+}{x}\right),`$ $`\left(\frac{\mathrm{\Delta }^{}}{x}\right)`$, $`\mathrm{\Delta }^+`$ and $`\mathrm{\Delta }^{}`$

$`\mathrm{\Delta }^+`$ $`=`$ $`M_{11}\mathrm{\Delta }^{}+M_{12}\left({\displaystyle \frac{\mathrm{\Delta }^{}}{x}}\right),`$ (32)
$`\left({\displaystyle \frac{\mathrm{\Delta }^+}{x}}\right)`$ $`=`$ $`M_{21}\mathrm{\Delta }^{}+M_{22}\left({\displaystyle \frac{\mathrm{\Delta }^{}}{x}}\right),`$ (33)
which was used by De Gennes in his general description of the Josephson junction close to $`T_c`$ . Following the approach , the critical current of the SFS junction is given directly by the expression

$$I_c=\frac{ie\pi N_s(0)D_s}{4T_cM_{12}}(\mathrm{\Delta }^{}\mathrm{\Delta }^+c.c),$$
(34)
where $`N_s(0)`$ is the electron density of states in S electrodes when they are in the normal state.

The calculation of the coefficient $`M_{12}`$ is straightforward but rather cumbersome for a general case. The corresponding formula is to intricate to be presented here. So, to discuss the essential physics, we consider below only the most interesting limits.

In the case of highly transparent interfaces, it is possible to obtain an analytical expression of the critical current taking into account the first order correction over $`h\tau `$ in the Usadel equation $`\left(\text{12}\right)`$

$$I_c=\frac{\pi \mathrm{\Delta }^2}{8eT_c}\frac{\sigma _n}{\xi _f}|\frac{k}{\mathrm{sinh}\left(\frac{kd}{2}\right)\mathrm{cosh}\left(\frac{kd}{2}\right)}+c.c.|$$
(35)
where the complex wave vector $`k=[(10.2h\tau )+i(1+0.2h\tau ]\sqrt{\frac{h}{D_f}}.`$ The dependences $`I_c(d)`$ for different parameter $`h\tau `$ are presented in Fig. $`2`$.

When the critical current goes through zero, the transition from $`0`$ to $`\pi `$-shift Josephson contact takes place. As expected, with the increase of the scattering time the period of oscillations becomes somewhat shorter, while its amplitude increases. This is consistent with a power low decrease of the critical current with F-layer thickness in the clean limit .

Further we concentrate on the influence of the S/F boundary transparency on the $`I_c`$ oscillations, and then we neglect the $`h\tau `$ correction in the diffusion coefficient. Starting with the case of completely transparent interfaces $`\gamma _{B1,2}0,`$ we retrieve the corresponding expression previously obtained in

$$I_c=\frac{\sigma _n}{\xi _f}\frac{\pi \mathrm{\Delta }^2}{eT_c}\left|\frac{\mathrm{sin}(x)\mathrm{cosh}(x)+\mathrm{cos}(x)\mathrm{sinh}(x)}{\mathrm{cos}(2x)\mathrm{cosh}(2x)}\right|,$$
(36)
where $`x=d/\xi _f.`$ In the case of low transparency of both interfaces $`\gamma _{B1,2}1,`$ the critical current is

$$I_c=\frac{1}{\gamma _{B1}\gamma _{B2}}\frac{\sigma _n}{\xi _f}\frac{\pi \mathrm{\Delta }^2}{2eT_c}\left|\frac{\mathrm{cos}(x)\mathrm{sinh}(x)\mathrm{sin}(x)\mathrm{cosh}(x)}{\mathrm{cos}(2x)\mathrm{cosh}(2x)}\right|.$$
(37)
If one interface (the left S/F interface) is transparent, $`\gamma _{B1}0,`$ the critical current is

$$I_c=\frac{\sigma _n}{\xi _f}\frac{\pi \mathrm{\Delta }^2}{eT_c}\left|\frac{\mathrm{cosh}(x)\mathrm{sin}(x)+\mathrm{cos}(x)\mathrm{sinh}(x)+2\gamma _{B2}\mathrm{cos}(x)\mathrm{cosh}(x)}{(12\gamma _{B2}^2)\mathrm{cos}(2x)(1+2\gamma _{B2}^2)\mathrm{cosh}(2x)2\gamma _{B2}\left[\mathrm{sin}(2x)+\mathrm{sinh}(2x)\right]}\right|$$
(38)
If one interface (the first S/F interface) has a low transparency, $`\gamma _{B1}1,`$ the critical current is

$$I_c=\frac{\sigma _n}{\gamma _{B1}\xi _f}\frac{\pi \mathrm{\Delta }^2}{eT_c}\left|\frac{\mathrm{cos}(x)\mathrm{cosh}(x)\gamma _{B2}\left[\mathrm{cosh}(x)\mathrm{sin}(x)\mathrm{cos}(x)\mathrm{sinh}(x)\right]}{(12\gamma _{B2}^2)\mathrm{cos}(2x)+(1+2\gamma _{B2}^2)\mathrm{cosh}(2x)2\gamma _{B2}\left[\mathrm{sin}(2x)\mathrm{sinh}(2x)\right]}\right|$$
(39)
The evolution of the $`I_c(d)`$ dependence with the decrease of the transparency of the second interface is presented in Fig. $`3`$. We observe that with the increase of $`\gamma _{B2}`$ the amplitude of oscillations decreases, as well as the thickness of the ferromagnetic layer corresponding to the first zero of the critical current.

In the experiment , the presence of Al2O3 barrier at one S/F boundary can be modeled by a low transparency interface $`\left(\gamma _{B1}1\right)`$ while the other boundary is quite transparent $`\gamma _{B2}0.`$ The resistance $`R_n`$ per unit area of the S/F/S junction in this experiment is dominated by the resistance of the tunnel barrier and thus $`R_n\frac{\gamma _{B1}\xi _f}{\sigma _n}`$. So we may expect that the experimental situation  must be described by the expression $`\left(\text{39}\right)`$ with $`\gamma _{B2}=0`$ i.e.

$$I_cR_n=\frac{\pi \mathrm{\Delta }^2}{eT_c}\left|\frac{\mathrm{cos}(x)\mathrm{cosh}(x)}{\mathrm{cos}(2x)+\mathrm{cosh}(2x)}\right|.$$
(40)
Note that this expression is quite different from that used in  to fit the experimental data. Unfortunately in  there is no information which could shed light not only on the reasons of the discrepancy between our formula but even at the origin of the used theoretical expression. The experimental results on the $`I_c(d)`$ oscillations presented in  have been obtained at low temperature. Nevertheless we think that the overall shape of the $`I_c(d)`$ curve is not very sensitive to the temperature and we compare the experimental data  with our theoretical expression $`\left(\text{40}\right)`$ in Fig. $`4.`$ We use as the fitting parameter $`\xi _f30`$ Å, while the experiment  provided a value of $`\xi _f`$ around $`35`$ Å in our notation. The value of the other fitting parameter $`\frac{\pi \mathrm{\Delta }^2}{eT_c}`$is $`110`$ $`\mu V`$; this parameter is quite difficult to estimate because of the uncertainty on the value of $`\mathrm{\Delta }`$ at the S/F interface in the geometry of experiment . Also following the analysis of the authors  presented in their previous publication, we have taken into account that the actual ferromagnetic thickness of PdNi layer is reduced by $`15`$ Å due to some interdiffusion at the S/F interface. The obtained description of the experimental data  is quite satisfactory.

In conclusion, we have presented the general theoretical description of the ferromagnetic $`\pi `$-junctions near superconducting transition temperature and proposed a simple way to take into account the finite elastic scattering time in Usadel equation. The obtained analytical expressions may be useful for the analysis of different experimental realizations of such junctions. Our analysis may be also easily generalized to the situation when the superconducting electrodes are fabricated from different materials.

We thank M. Aprili and V. Ryazanov for stimulating discussions. This work was supported by the ESF “vortex” Programme and the ACI “supra-nanométrique”.

FIG. $`1`$ Geometry of the considered S/F/S system. The thickness of the ferromagnetic layer is $`d`$. The transparency of the left S/F interface is characterized by the coefficient $`\gamma _{B1}`$ and the transparency of the right F/S interface is characterized by $`\gamma _{B2}.`$

FIG. $`2`$ Critical current of the S/F/S junction as a function of the thickness of the ferromagnetic layer normalized by $`\xi _f`$ for different scattering time. Both S/F interfaces are completely transparent. The parameter $`h\tau `$ is equal to $`0`$, $`0.25`$ and $`0.5`$.

FIG. $`3`$ Critical current of the S/F/S junction as a function of the thickness of the ferromagnetic layer normalized by $`\xi _f`$. The first S/F interface has a low transparency $`\left(\gamma _{B1}1\right)`$. The parameter $`\gamma _{B2}`$ characterizing the transparency of the second interface is chosen as $`0.2,`$ $`1.5,`$ $`3.0,`$ and $`10.`$

FIG. $`4`$ The experimental points correspond to the measurement of the critical current, done by Kontos et al , vs the PdNi layer thickness. The theoretical curve is the best fit obtained by using formula $`\left(\text{40}\right)`$. The fitting parameters are $`\xi _f30`$ Å and $`\frac{\pi \mathrm{\Delta }^2}{eT_c}110`$ $`\mu V.`$


text:
# Untitled Document

Fractional Charge in Quantum Hall effect Keshav N. Shrivastava School of Physics, University of Hyderabad,
Hyderabad 500046, India

In 1976 Jackiw and Rebbi found 1/2 of a fermion number by using Dirac equation. Schrieffer, in several proposals made an effort to suggest that there is a fractional charge. The calculations of Peierls distortion, Berry’s phase and classical action were presented to accomodate the fractional charge in non-relativistic theory. Laughlin used the antisymmetry to define the charge density per unit area in a two dimensional system. In order to elliminate the area, Laughlin introduced the incompressibility which fixed the area and the odd number which determines the antisymmetry of the electron wave function, gave the charge. The antisymmetry relies on the odd number which can be equal to 3 so that the charge became 1/3. We have used the orbital angular momentum and the spin of the electron to define the effective charge through the Bohr magneton which already has the charge of the electron. When both signs of the spin are used, $`j=l\pm s`$, we generate a series of charges which are in agreement with those measured by Stormer in the quantum Hall effect experiments.

Corresponding author: keshav@mailaps.org
Fax: +91-40-2301 0145.Phone: 2301 0811.

1. Introduction

In 1976, Jackiw and Rebbi, while considering the expansion of the Fermi quantum fields in terms of eigen functions of the Dirac equation, introduced zero-energy solution in addition to positive and negative energy solutions, in 1+1 dimensions. They obtained conserved Fermi number current but the fermion number is found to be +1/2 for the soliton and -1/2 for the antisoliton. Later, Su and Schrieffer decided to introduce 1/2 the charge of an electron in systems which undergo Peierl’s distortion. Considerable effort was put to find the fractional charge in polyacetelene, $`CH_x`$ chain. It was thought that $`CH_x`$ chain will exhibit fractional charge but the ferromagnetic state of such a fractionally charged state is not known. In 1983, Laughlin wrote a wave function on the basis of antisymmetry from which it was suggested that the quasiparticles can have 1 or 1/3 charge. Since antisymmetry requires that the exponent should be an odd number in order to generate a negative sign in the wave function with respect to interchange of any pair of particles. When this odd number is +1, interchanging a pair of particles changes the wave function $`\psi `$ to $`\psi `$. For example, $`z_1z_2`$ changes sign when 1 and 2 are interchanged. If the exponent is 3, we have $`\psi (z_1z_2)^3`$, which upon exchanging 1 and 2 also changes $`\psi `$ to $`\psi `$. Laughlin defined the charge of the quasiparticles in such a way that the odd exponent 1 or 3 became the charge of the quasiparticles 1 or 1/3. Thus, a wave function was generated in which excitations will have a charge of 1/3. It was argued that this determines the wave function of quasiparticles, which cause steps in the quantum Hall effect experiments in semiconductors. Since the exponent must be an odd integer, the fractional charges produced are 1/3, 1/5, 1/7, etc. There was considerable belief that such charges are indeed found in the experimental data. If the charge became 1/3, what happened to the remaining charge of 2/3 is not discussed. It is not clear whether there is a phase transition from e to e/3. In the Jackiw and Rebbi’s theory the missing charge is carried by the soliton but in the Laughlin’s theory one is concerned only with the wave function of the fractional charge.

While several authors were interested in seeing, if there is a spin-charge decoupling, we have found in 1986 that spin and charge are in fact coupled. This means that the knowledge of spin is needed to determine the value of the charge. The charge can be written in terms of spin and orbital angular momenta quantum numbers, $`l`$ and $`s`$,as,

$$\frac{e_{eff}}{e}=\frac{l+\frac{1}{2}\pm s}{2l+1}$$
(1)
It has its origin in the factors which multiply the Bohr magneton, $`\mu _B=e\mathrm{}/2mc`$ so that the charge $`e`$ can be replaced by $`e_{eff}`$. Usually, there is only one $`g`$ value associated with the Lande’s formula which uses only one sign in L+S but here we are using half the value with both signs, $`j=l\pm s`$. When the charge $`e_{eff}=0`$, we see that there is a particle with spin $`\pm 1/2`$. The zero charge also corresponds to zero energy, so that at this point there is a charge-density wave. Thus we predict a particle which is electrically neutral but has a spin. Not only that, but there is an effort to identify these objects in the experimental data of quantum Hall effect and it is claimed that such quasiparticles really exist.

In this paper, we describe the suggestions of Schrieffer about the fractional charges, present a brief review of Laughlin’s wave function and our own work on the fractional charges.

2. Schrieffer’s suggestions

(a) Peierls distortion.

It has been suggested by Su and Schrieffer that in the case of dimerization, if the electrons were spinless, the charge of a kink will be $`+\frac{1}{2}|e|`$. In the case of spin, charged kinks are spinless and neutral kinks have spin $`1/2`$ in contrast to $`\pm e`$ charge for electrons and holes in semiconductors and spin $`1/2`$, i.e.,

$$charge=\pm \frac{1}{2}e,spin=0;$$
(2)
and

$$charge=0,spin=1/2;$$
(3)
and the usual,

$$charge=\pm e,spin=1/2;$$
(4)
excitations are possible. In the one-third filled band case the charge deficit is $`(1/3)e`$ or $`(2/3)e`$ and one is left with a fractional charge,

$$Q=\pm \frac{1}{3}|e|or\frac{2}{3}|e|,$$
(5)
but there is no way of knowing, what determines whether there is dimerization or trimerization. For commensurability $`n`$, the kinks have charge which is an integer multiple of $`\pm e/n`$.

Consider three types of distortions with sites numbered from -1, 0, 1, 2, …,8. The dimerizations occur as (-1,0), (2,3), (5,6) which means that atoms on sites -1 and 0 are paired, those on 2 and 3 are paired and (5,6) are paired. This is called A dimerization. Similarly (0,1), (3,4), (6,7) is called B and (1,2), (4,5), (7,8) is called C. There are two types of kinks, type I is A to B, B to C or C to A and type II are A to C, C to B or B to A. The primitive charge of type-I kink is $`+(2/3)|e|`$ and that of type-II kink is $`(2/3)|e|`$. In the case of type-I kinks, there is a localized state, $`\phi _l`$ and symmetrically a state $`\phi _u`$ in the lower half of the upper gap. If $`\phi _l`$ has a charge of $`(2/3)|e|`$ with all spins paired, the kink spin is zero. If $`\phi _l`$ is slightly occupied, the kink has charge $`(1/3)|e|`$ and spin 1/2. If $`\phi _l`$ is doubly occupied, the charge is $`(4/3)|e|`$ and the spin is zero. For type-II kinks, the localized states $`\phi _l`$ and $`\phi _u`$ are below and above the gap centers, respectively. If $`\phi _l`$ is unoccupied, the kink has charge $`+(4/3)|e|`$ and spin 0. For $`\phi _l`$, singly occupied, the charge is $`(1/3)|e|`$ and spin is (1/2), while for $`\phi _l`$ doubly occupied, the charge is $`(2/3)|e|`$ and spin is zero. In other words, fractional charges, with or without spin occur.

In the above discussion, we have not considered the spin=0, charge=0; spin=0, charge $`\pm `$ 1 or spin = 1/2, charge = $`\pm `$ 1/2. The spin = zero, charge =0 is perfectly correct because it is a neutral boson; spin=0, charge =$`\pm `$ 1 is a charged particle without spin which is also very reasonable and spin =1/2, charge = $`\pm `$ 1/2 has unusual symmetry. However, due to Goldstone theorem, only spin=0 quasiparticle should emerge in a Peierls distortion. Therefore, the fractional charge which may not be a boson violates the Goldstone theorem.

(b) Berry’s phase.

Arovas et al calculate the change of phase $`\gamma `$ of a state having a quasihole localized at $`z_o`$. As $`z_o`$ adiabatically moves around a circle of radius R enclosing flux $`\varphi `$, the accumulated phase factor is,

$$\frac{e}{\mathrm{}c}𝐀.d𝐥=2\pi (\frac{e}{e})(\frac{\varphi }{\varphi _o}).$$
(6)
This is the phase gained by the quasiparticle of charge e\* in going around a loop. It may be noted that whether e\* is changed or A is changed, will be indistinguishable since only the product e\*A occurs in the phase factor. It will be desirable to obtain a link with the wave functions suggested by Laughlin for the quasiparticles of fractional charge. Laughlin’s calculation is concerned with the charge per unit area and not only the quasiparticle charge so that e\* should be replaced by e\*/m$`a_o^2`$, where m is an odd integer and $`a_o^2`$=$`\mathrm{}c/eB`$, which involves the magnetic induction, B. The charge $`e`$ in absence of a field is linked to the Bohr radius, $`a_B`$. The corrected form of the phase factor must have $`ma_o^2`$. Since, the phase must be dimensionless, the corrected form of the phase should be,

$$\gamma =2\pi (\frac{e}{e})(\frac{a_B^2}{ma_o^2})(\frac{\varphi }{\varphi _o}).$$
(7)
This means that the phase factor is linearly dependent on the magnetic field. At the field of a few Tesla, $`a_o10^4cm`$,while $`a_B10^8cm`$ so that $`(a_B/a_o)^210^8`$. Therefore, the phase factor $`\gamma `$ is a negligibly small number. If $`\omega `$ is the frequency and $`\pi /a_o`$ is the wave vector, then the velocity is $`\omega a_o/\pi `$ and $`\gamma `$ is not independent of the velocity. The charge of the quasiparticle is $`\pm e=e/m`$ so that it is possible to obtain a fractional charge,

$$Q=\pm \nu e.$$
(8)
as a possibility.

(c) Classical action.

Kivelson et al have calculated the classical action and found that when the real part is less than a critical value, there is a phase transition. The fractional values of the charge such as $`e/3`$ satisfy the condition of a phase transition so that a phase transition is possible for a fractional charge. Here also we have found that whether $`e`$ is being changed to $`e`$ or $`l_o`$ is being changed to $`l_o/\nu ^2`$, are indistinguishable.

3. Laughlin’s theory

Laughlin utilized the antisymmetry property of the wave function of fermions to define a fractional charge. Let us consider a wave function of the form,

$$\psi =(z_1z_2)^mf_s$$
(9)
where $`f_s`$ is a function which is symmetric with respect to interchange of 1 and 2, i.e., there is no change in $`f_s`$ when the particles 1 and 2 are interchanged, $`f_s`$(1,2)=$`f_s`$(2,1). This property is not associated with the anticommutators for fermions. On the other hand, when m=1, the factor $`z_1z_2`$ changes sign when we interchange 1 and 2. Therefore, the factor $`z_1z_2`$ is associated with the antisymmetry. Usually, it is enough to have m=1 but it is clear that the antisymmetry is obtained as long as $`m`$ is an odd integer. For example, for $`m`$=3 also the factor containing $`z_1z_2`$ will change sign with respect to interchange of 1 and 2. The factor $`(z_1z_2)^m`$ is therefore associated with the antisymmetry of the wave function with respect to interchange of 1 and 2 as long as $`m`$ is an odd integer. Laughlin has found that the fluid density is,

$$\rho _m=\frac{1}{m2\pi a_o^2}$$
(10)
so that the charge density is $`e/[m2\pi a_o^2]`$. Laughlin interpreted that the charge has become $`e/m`$ or e/1, e/3, e/5, etc. For this interpretation, it is possible to explain the experimental data on noise measurements as well as plateaus in the quantum Hall effect. An equally valid interpretation is that charge does not fractionalize but $`a_o^2`$ does. That means, instead of $`e/m`$ and $`2\pi a_o^2`$, we select $`e`$ and $`2\pi (ma_o^2)`$ in the above expression for density. The antisymmetry property or the exactness of Laughlin’s calculation remain unaffected by this change in the interpretation. The square of the magnetic length is given by,

$$a_o^2=\frac{\mathrm{}c}{eH_o}.$$
(11)
Therefore, $`ma_o^2`$ has the magnetic field as the multiplier. Therefore, the quantity which appears in the charge density becomes $`ma_o^2H_o`$. Therefore, $`m`$ need not be attached to the charge. The energy of the electrons in a magnetic field is given by $`g\mu _BH.S`$ where charge $`e`$ is included in $`\mu _B`$. Therefore, four distinct quantities are involved. If any one of these is changed, the blame can easily go to any of the other three.

4. Theory.

We\[9-15\] use a vector model of $`𝐥`$ and $`𝐬`$ to define $`𝐣`$. In the Lande’s formula only $`𝐉=𝐋+𝐒`$ is used so that there is only one value but we use $`𝐣=𝐥\pm 𝐬`$ so that there are two values,$`g_\pm `$ which is different from Lande’s. When we define the Bohr magneton, the effective charge of the electron becomes,

$$e_{eff}=\frac{l+(\frac{1}{2})\pm s}{2l+1}e$$
(12)
For various values of $`l`$ = 0, 1, 2, 3, … $`\mathrm{}`$, this expression gives two series of charges, one for + sign and the other for minus sign. This predicted series is exactly the same as that measured by Stormer in the quantum Hall effect. The magnetic moment of the quasiparticles is given by,

$$\mu _{eff}=\frac{e_{eff}\mathrm{}}{2mc}.$$
(13)
Therefore, when we go from the classical Hall effect to the quantum Hall effect, the charge of the quasiparticles changes from $`e`$ to $`e_{eff}`$. A large amount of data on the quantum Hall effect has been examined by us and our theory is found to be correct in all the cases.

5. Conclusions.

In the case of an incompressible system $`a_o`$ has been left out so that the effect of $`m`$ is to change the charge to $`e/m`$ where $`m`$ is an odd integer associated with the antisymmetry of the wave function. On the other hand, we can keep $`e`$ unchanged and throw the effect of $`m`$ on $`a_o^2`$ so that the product $`e/ma_o^2`$ is unchanged. Since $`a_o`$ depends on the magnetic field, $`e/(ma_o^2)=e^2H_o/(m\mathrm{}c)`$, changing $`e`$ is equivalent to changing $`H_o`$. Hence, the change in $`H_o`$ depends on $`l`$ and $`s`$ and we obtain spin dependent flux quantization. We also obtain a quasiparticle which has zero charge and spin =1/2.

6. References

1. R. Jackiw and C. Rebbi,Phys. Rev. D13, 3398(1976).
2. W. P. Su and J. R. Schrieffer, Phys. Rev. Lett. 46,738(1981).
3. D. Arovas, J. R. Schrieffer and F. Wilczek, Phys. Rev. Lett. 53, 722 (1984).
4. R. B. Laughlin, Phys. Rev. Lett. 50, 1395 (1983).
5. K. N. Shrivastava, cond-mat/0211351.
6. S. Kivelson, C. Kallin, D. P. Arovas and J. R. Schrieffer, Phys. Rev. Lett. 56, 873 (1986).
7. K. N. Shrivastava, cond-mat/0211621.
8. K. N. Shrivastava, cond-mat/0210238.
9. K. N. Shrivastava, Natl. Acad. Sci. Lett. (India) 9,145(1986).
10. K. N. Shrivastava, Phys. Lett. A113, 435 (1986);A115,459(1986)(E).
11. K. N. Shrivastava, Mod. Phys. Lett. B 13,1087(1999).
12. K. N. Shrivastava, Mod. Phys. Lett. B 14, 1009(2000).
13. K. N. Shrivastava, p.235 in, Frontiers of Fundamental Physics 4, edited by B. G. Sidharth and M. V. Altaisky, Kluwer Academic, New York 2001.
14. K. N. Shrivastava, Superconductivity: Elementary Topics, World Scientific, New Jersey 2000.
15. K.N. Shrivastava, Introduction to quantum Hall effect,
Nova Science Pub. Inc., N. Y. (2002).
16. H. L. Stormer, Rev. Mod. Phys. 71, 875 (1999).

Note: Ref.15 is available from:
Nova Science Publishers, Inc.,
400 Oser Avenue, Suite 1600,
Hauppauge, N. Y.. 11788-3619,
Tel.(631)-231-7269, Fax: (631)-231-8175,
ISBN 1-59033-419-1 US$`\$69`$.
E-mail: novascience@Earthlink.net


text:
# Dispersion of imbibition fronts

11institutetext:  Institut Français du Pétrole, 92852 Rueil-Malmaison Cedex, France
Laboratorio de Física Estadística de Sistemas Desordenados, Centro de Física, IVIC, Apartado 21827, Caracas 1020-A, Venezuela
Center for Bioinformatics, Saarland University, P.O. Box 151150, 66041 Saarbrücken, Germany

Flows in porous media Diffusion; interface formation Wetting

## Abstract

We have studied the dispersive behaviour of imbibition fronts in a porous medium by X-ray tomography. Injection velocities were varied and the porous medium was initially prewetted or not. At low velocity in the prewetted medium, the imbibition profiles are found to be distinctly hyperdispersive. The profiles are anomalously extended when compared to tracer fronts exhibiting conventional (Gaussian) dispersion. We observe a strong velocity dependence of the exponent characterizing the divergence of the dispersion coefficient for low wetting-fluid saturation. Hyperdispersion is absent at high imbibition velocities or when the medium is not prewetted.

Dispersion is a process by which an initial distribution of fluid spreads out in a porous medium under the effects of the disordered nature of the porous structure and the complex flow events occurring in the pores. This process has practical relevance in the fields of hydrology (water infiltration, propagation of pollutants in soils ), petroleum engineering (oil recovery by water injection ), and biotechnology (gel and chromatographic analysis ).

A dispersion process can be described on the macroscopic (Darcy) scale, i.e., for an average fluid content (or saturation) $`S`$ defined over a “representative elementary volume” , by the equation

$$\varphi _tS+u_xS=_x(D_xS)$$
(1)
in which $`u`$ is the average velocity (or drift), $`x`$ the position along the flow direction, and $`D`$ is the dispersion coefficient. The simplest dispersion process, referred to as hydrodynamic dispersion, concerns two fluids of identical nature in which one contains “tagged” or tracer molecules. Then, $`S`$ is the average tracer concentration, and dispersion arises primarily from the variation in local fluid velocities (molecular diffusion is negligible for practical values of $`u`$). For homogeneous and non-consolidated porous media, the dispersion coefficient is given by $`D=u\mathrm{}`$, where $`\mathrm{}`$ is a length of the order of the typical pore size. Dispersion is Gaussian (i.e., $`D`$ is independent of $`S`$), and $`S`$ is a simple function of the reduced variable $`(xut/\varphi )/\sqrt{t}`$ . For many “real” porous media, however, one observes a non-Gaussian hydrodynamic dispersion: the physical mechanisms governing such anomalous behaviour are being actively investigated .

Here we are concerned with the dispersion behaviour found in imbibition: a front of wetting fluid is injected through one face of a horizontally placed porous medium. The medium is initially filled with a nonwetting fluid. Therefore, a front between wetting and nonwetting fluid forms which spreads as it advances into the porous medium, as illustrated in Figure 1.

In this case, dispersion is primarily caused by capillary forces which, depending on pore size and structure, accelerate or slow down the advance of the wetting fluid: it is usually referred to as capillary dispersion . Imbibition is governed by an equation similar to eq.(1) in which the convective term is replaced by $`u_xF(S)`$, where $`F(S)`$ is the (generally nonlinear) fractional flow function from which dispersion-free solutions are constructed. They are travelling waves or shocks, depending on the functional form of $`F`$ and the initial and boundary conditions . These solutions are good approximations of the true solutions when convective (viscous) forces dominate capillary forces, i.e., for high capillary numbers.

Anomalous (non-Gaussian) dispersion behaviour has been suggested theoretically to occur  and observed experimentally  in imbibition at low or vanishing injection velocity (‘passive’ imbibition). To characterize the evolution of the imbibition fronts a power-law divergence of the dispersion coefficient at low wetting fluid saturation is used. In the proposed models  the (velocity-independent) exponent is related to the pore structure (through a fractal dimension) and to the range of the substrate-fluid interactions.

In this letter, we present and discuss a series of imbibition experiments carried out at different imposed velocities and initial wetting states of the porous medium. The experiments consist in monitoring the progression of the front of the wetting phase (water) displacing a nonwetting phase (oil) by means of an X-Ray CT-Scanner (CT = computerized tomography). From the measured water saturation profiles we are able to characterize the dispersion behaviour in a precise manner, especially at low injection velocity or capillary number, as we will show. We note that for our system, we have no a priori knowledge of either the functional dependence of $`F(S)`$ nor of the dispersion coefficient $`D`$ and its possible dependence on wetting fluid saturation .

Our experiments were performed in a porous medium made of quartzitic grains with radii $`100\mu `$m, packed in a cylindrical glass column with a diameter of $`3`$ cm, and length $`L=40`$ cm. The porosity (or void percentage) and absolute permeability were $`\varphi 35\%`$, and $`k10\mu `$m2 (or Darcy), respectively. Homogeneity was checked from X-Ray scanner measurements of the dry porous medium. The X-Ray CT-Scanner device used in this study was a fourth generation medical model FX/i, General Electric that allow very rapid measurements: less than 1 minute suffices to measure about 40 images of cross-sections of the porous medium, each cross-section being placed every cm apart. As the displaced fluid we used n-decane with 99$`\%+`$ purity, and the (distilled) water contained a dopant (30 g/l potassium iodide) to enhance its X-Ray contrast with the oil. The displacing fluid (water) is slightly more viscous than the displaced fluid (oil), 1.03 vs. 0.89 $`mPas`$, and their interfacial tension at the temperature of the experiments (23 C) is $`\sigma 53`$ dyn/cm. The capillary number, i.e. the ratio of viscous to capillary forces $`Ca=\mu u/\sigma =\mu Q/\gamma A`$ spans more than two orders of magnitude, from around 10-8 to 10-6. Here, $`Q`$ is the imposed injection rate, and $`A`$ the cross-sectional area of the porous medium. The mobility ratio between the two fluid phases is favorable and one therefore expects the development of a front which is stable with respect to viscous fingering.

The procedure for our experiments in a prewetted porous medium is as follows. First the dry medium is fully saturated with water. Then, it is flooded at a high rate with an excess of oil, until no water leaves the medium. Some water cannot be displaced and remains trapped in the medium. The profile $`S_{wi}(x)`$ of this irreducible water measured with the CT-Scanner amounts to 25$`\%`$ all across the medium, except very near the outlet where it reaches 45$`\%`$. Such an accumulation of water is nothing but the capillary end effect. Then the experiment starts by injecting the water at an imposed flow rate $`Q`$. Water saturation profiles $`S_w(x,t)`$ are measured at various times $`t`$ following the start of the injection. The $`S_w(x,t)`$ are calculated from the measured data by linear interpolation with previous CT-Scan data of the porous medium fully saturated with water or oil.

For an imbibition experiment conducted at a different rate, the porous medium is again first flooded with oil (as described above) to reach the state of irreducible water saturation $`S_{wi}(x)25\%`$. Imbibition experiments in a prewetted medium were conducted at the flow rates $`Q=`$ 1, 3, 10, 30, 100, and 300 ml/hr, corresponding to superficial velocities $`u/\varphi `$ (these velocities would correspond to those of sharp advancing fronts) which equal $`0.14,0.42,1.4,4.2,14`$, and $`42`$ cm/hr ($`Ca`$ ranging from $`0.8`$ $`10^8`$ to $`2.4`$ $`10^6`$).

We now turn to the description of our findings. The saturation profiles corresponding to the two lowest velocities in a prewetted porous medium are shown in Figs. 2. At low injection rates (low $`Ca`$) the profiles increase almost uniformly throughout the entire porous medium. This behaviour qualitatively indicates the presence of a hyperdispersive imbibition process governed by a dispersion coefficient which diverges at low saturations as a power law . This behaviour is analyzed quantitatively below.

We first assess the role of prewetting on the observed hyperdispersive behaviour, which we infer from a comparison to a reference experiment in a non-prewetted porous medium. In this case, the medium was first fully saturated with oil and then the water was injected at the imposed velocity. The resulting profiles are shown in Fig. 3. The saturation front is well-defined and localized, and the signature of hyperdispersion displayed by the fronts in Fig. 2 is clearly absent.

By comparison, Figs. 4 a), b) show the saturation profiles in the medium in the prewetted medium at high velocity. With the increased velocity (i.e., larger $`Ca`$) the imbibition fronts sharpen and become very straight at the two highest rates ($`Q=100`$ and 300 ml/hr). Such a behaviour might be indicative of a hypodispersive behaviour, i.e. $`D0`$ for $`S_w0`$ .

In order to further characterize the behaviour of these fronts we have conducted tracer experiments used as a probe for Gaussian (hydrodynamic) dispersion. In those tracer experiments, the water doped with the tracer molecules (30 g/l potassium iodide) was injected at a fixed rate which displaced the pure water initially fully saturating the porous medium. Then the tracer concentration profiles were measured at various successive times $`t`$ with the X-Ray CT-Scanner. We checked that the tracer profiles could be superimposed when expressed in terms of the reduced variable $`(xut/\varphi )/\sqrt{t}`$ . In Fig. 4(a) and (b), the data sets of the tracer fronts are superimposed with the highest flow rate imbibition fronts. Unfortunately we have too few data to quantify the capillary dispersion coefficient at low $`S_w`$ to unambiguously decide on the presence of hypodispersion in the fronts. On a longer distance (more than the 40 cm of the porous medium available to us) we expect the profiles to turn into shocks .

We now turn to the quantitative analysis of the saturation profiles at low flow velocities. In hyperdispersion, the dispersion coefficient diverges at low wetting phase saturation according to a power law . Early experiments by Bacri et al.  had previously shown evidence for this effect in a prewetted porous medium. In their case, the data were analyzed assuming $`DS_w^1`$, but the acoustic technique used allowed saturation measurements only for very few positions along the porous medium.

Based on our data, we indeed find that the saturation profiles can be interpreted in terms of a power law divergence of the dispersion coefficient, however of the form

$$D(S_w)(S_wS_{wi})^n,n<0$$
(2)
in which case eq.(1) has a self-similar scaling solution for the saturation profile for small $`u`$ . Consequently, the saturation profiles superimpose when scaled by $`t^{1/(n+2)}`$ if positions are rescaled by $`t^{(n+1)/(n+2)}`$. By using this scaling solution, we can collapse our low flow rate saturation data onto a single master curve (Fig. 5).

However, the exponent $`n`$ allowing the best collapse exhibits a strong dependence with the imposed flow rate: $`n1.1`$ for $`Q=3`$ ml/hr, and $`n0.5`$ for $`Q=10`$ ml/hr. For the lowest rate $`Q=1`$ ml/hr, where a “giant hyperdispersivity” is found characterized by an almost homogeneous profile spanning the whole system (see Fig. 2 a), we could not determine a scaling exponent.

These results have not been anticipated theoretically. In Refs.  it was argued that a divergence of $`D`$ should be related in a simple scaling fashion to the fluid-substrate molecular interactions that determine the thickness of wetting films, or to surface geometry. Our porous medium does not possess a complex (fractal) pore geometry ; in addition, since the wetting fluid occupies about 25$`\%`$ of the pore space, the thickness of the wetting layers is much larger than the range (several tens of nanometers) of the fluid-substrate interactions. In our experiments it is rather capillarity that controls the distribution of the fluid in the pores space. Most of the wetting phase is in the “bulk” form of menisci and “pendular” structures, sometimes referred to as layers , occupying corners and substrate concavities. These wetting layers provide a continuous and conductive pathway across the porous medium, and presumably swell rather uniformly across the porous medium at low injection velocities.

A challenge for future studies is to understand how the flow events on the pore scale give rise to the observed hyperdispersive behaviour on the macroscopic (Darcy) scale. These events are complex, as they include cooperative pore invasion, flow in layers along the substrate micro-roughness, and collapse between wetting layers. A very useful tool for modeling these events is provided by pore-network models that have shed some light on various patterns and macroscopic flow parameters observed in imbibition under variation of capillary number, contact angle and wetting phase saturation . However, so far these models do not account for a variation in the conductance of the wetting layers (i.e., their swelling) which we believe is needed to reproduce the strong hyperdispersive effect we observe. Some progress in pore-network modeling is to be expected from the incorporation of proper swelling rules, e.g. for layers in corners .

To conclude, we have demonstrated experimentally, for the first time, that the presence of wetting layers is needed for hyperdispersion to occur in imbibition in a porous medium. For low wetting-phase saturation the dispersion coefficient behaves like a power-law with a velocity-dependent effective exponent. New theoretical advances are needed to account for these results which we expect can serve as a benchmark for future pore-network modeling studies.

###### Acknowledgements.

We are grateful to C. Fichen and C. Schlitter for their assistance in the CT-Scan measurements, and to O. Vizika and B. Noetinger for discussions. Y.M. was supported by the Post-Graduate Cooperation Program (PCP) France-Venezuela. R.B. was supported by DFG under Schwerpunktprogramm “Wetting and Structure Formation at Surfaces” (Bl-256) and the Flemish Government under Grant (VIS/97/01). We further acknowledge support under TOURNESOL grant T2000.013.