text:
# Domain walls of relative phase in two-component Bose-Einstein condensates

## I Introduction

One of the long-sought goals in low-temperature physics is the creation of two interpenetrating superfluids. Early efforts were directed at mixtures of helium isotopes. More recently, following the experiments with Bose-Einstein condensates (BEC) of atomic gases , considerable efforts have been made to create systems where two species of atoms condense simultaneously. This goal was partially achieved for two different hyperfine spin states of 87Rb, which were condensed in the same trap by the technique of sympathetic cooling . Later the dynamics of the complex relative motion of the condensates has been studied . The possibility of the measurement of the relative phase between the two condensates has also been demonstrated . In these experiments the two condensates have a substantial overlap, although they do not completely interpenetrate each other in the stationary state. A similar state called “spinor condensate” has been observed for sodium gas .

Theoretical investigation of two-component Bose systems has started many decades ago with the construction of the phenomenological hydrodynamic equations in the spirit of the Landau-Khalatnikov two-fluid model for the one-component BEC . Later, this construction has been put onto a microscopic basis . Recent experiments with alkali atoms have revived the interest in the theory of such systems. Hartree-Fock theory has been successfully tested on the two-component 87Rb system . The stability , ground-state properties  and collective excitations  have been studied theoretically by using the Gross-Pitaevskii equations.

Many properties of two-component, or binary, BEC can be understood from symmetry arguments. Compared to one-component Bose superfluids, two-component systems have more interesting pattern of symmetry and symmetry breaking. Condensation in binary Bose systems corresponds to the spontaneous breaking of two (instead of one) global U(1) symmetries. These symmetries are related, by Noether’s theorem, to the separate conservation of the number of atoms of each of the two species. The quantum state of the binary Bose system, therefore, is characterized by two phases of the two condensates. Correspondingly, the physics of binary BEC is also richer than of usual one-component systems.

The effects of a symmetry are often best exposed by violating the symmetry explicitly in a controlled fashion. A very interesting feature, specific to systems consisting of atoms of the same isotope in different spin states, is that it is possible to couple two condensates by a driving electromagnetic field tuned to the transition frequency. In this case atoms can be interconverted between the two spin states and the numbers of atoms of each species are not conserved separately anymore; only the total number of atoms is constant. This implies that, in the presence of the coupling drive, only one U(1) symmetry remains exact, the other one is explicitly violated. The preserved U(1) symmetry obviously comes from the conservation of the total number of atoms, and corresponds to changing the phases of the two condensates by the same amount (i.e., leaving the relative phase unchanged). The violated U(1) corresponds to changing the relative phase between the two condensate. The presence of the coupling drive lifts the degeneracy of the ground state with respect to the relative phase.

In this work, we show that a sufficiently small violation of the U(1) symmetry corresponding to the relative phase leads to the existence of a nontrivial static configuration — a domain wall inside which the relative phase changes by $`2\pi `$. This configuration is a local minimum of the energy. However, the domain wall is not topologically stable and can “unwind” itself. To unwind, however, the system must overcome an energy barrier. Thanks to this fact, the rate of the spontaneous decay of the domain wall is exponentially suppressed.

Our paper is organized as follows. Section II introduces the field-theoretical description of binary BEC. In Sec. III we describe the domain wall configuration, whose physical interpretation is given in Sec. IV. Section V deals with the boundary of finite domain walls and the related phenomenon of “vortex confinement”. Section VI contains concluding remarks. In Appendix A we find the metastability condition for the domain wall in the particular case when the densities of the two components are equal, and in Appendix B two different mechanisms for the decay of the domain wall, operating at different temperature regimes, are considered.

## II The Lagrangian, its symmetries and normal modes

In this Section, we use field theory to describe general properties of binary BEC. Our goal is to introduce notations and the formalism to lay the ground for the discussion of the domain walls in the next Section.

A binary dilute Bose system is described by a quantum field theory of two complex scalar fields $`\psi _1`$ and $`\psi _2`$. These fields have the meaning of the wave functions of the two condensates. The dynamics of these fields is governed by the following Lagrangian,

$$L=i\mathrm{}(\psi _1^{}_t\psi _1+\psi _2^{}_t\psi _2)H(\psi _1,\psi _2),$$
(1)
where the Hamiltonian $`H(\psi _1,\psi _2)`$ has the form

$`H(\psi _1,\psi _2)`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{2m}}(|\psi _1|^2+|\psi _2|^2)\mu _1\psi _1^{}\psi _1\mu _2\psi _2^{}\psi _2`$ (3)
$`+{\displaystyle \frac{g_{11}}{2}}|\psi _1|^4+{\displaystyle \frac{g_{22}}{2}}|\psi _2|^4+g_{12}|\psi _1|^2|\psi _2|^2{\displaystyle \frac{1}{2}}\mathrm{}\mathrm{\Omega }(\psi _1^{}\psi _2+\psi _2^{}\psi _1).`$
In Eq. (1) $`\mu _{1,2}`$ are the chemical potentials of the two species,***In real experiments $`\mu _1`$ and $`\mu _2`$ are functions of coordinates. We assume here that the trapping potentials are sufficiently wide so that these chemical potentials can be put to constants. $`g_{ij}`$ is the scattering amplitude, in the zero momentum limit, between an atom of the $`i`$-th species and that of the $`j`$-th species, and are proportional to the scattering lengths $`a_{ij}`$,

$$g_{ij}=4\pi \mathrm{}^2a_{ij}/m;$$
(4)
and $`\mathrm{\Omega }`$ is the Rabi frequency arising from the coupling drive.

By varying the action $`S=𝑑t𝑑𝒙L`$ with respect to $`\psi _{1,2}`$, the familiar Gross-Pitaevskii equations are directly obtained:

$`i\mathrm{}_t\psi _1`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{2m}}^2\psi _1\mu _1\psi _1+(g_{11}|\psi _1|^2+g_{12}|\psi _2|^2)\psi _1{\displaystyle \frac{1}{2}}\mathrm{}\mathrm{\Omega }\psi _2;`$ (5)
$`i\mathrm{}_t\psi _2`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{2m}}^2\psi _2\mu _2\psi _2+(g_{12}|\psi _1|^2+g_{22}|\psi _2|^2)\psi _2{\displaystyle \frac{1}{2}}\mathrm{}\mathrm{\Omega }\psi _1.`$ (6)
Let us start by finding the ground state when the coupling drive is off, $`\mathrm{\Omega }=0`$. In the superfluid ground state, both $`\psi _1`$ and $`\psi _2`$ have nonzero expectation values. These can be found by minimizing the potential energy part in Eq. (1) with respect to $`\psi _1`$ and $`\psi _2`$. This minimization procedure gives the equations determining the densities $`n_1=|\psi _1|^2`$, $`n_2=|\psi _2|^2`$ in terms of the chemical potentials $`\mu _1`$ and $`\mu _2`$,

$`g_{11}n_1+g_{12}n_2`$ $`=`$ $`\mu _1,`$ (7)
$`g_{12}n_1+g_{22}n_2`$ $`=`$ $`\mu _2.`$ (8)
More conveniently, one could view Eq. (8) as the equations fixing the chemical potentials for given values of the densities.

Strictly speaking, Eqs. (8) only correspond to an extremum of the potential energy. For it to be a local minimum, the quadratic form $`g_{11}n_1^2+g_{22}n_2^2+2g_{12}n_1n_2`$ needs to be positive definite:

$$g_{11}g_{22}g_{12}^2>0.$$
(9)
In fact, Eq. (9) is the condition for the mixture of the two Bose superfluids to be thermodynamically stable against segregation . In this paper we shall assume that (9) is satisfied.

In principle, the constants $`g_{11}`$, $`g_{12}`$ and $`g_{22}`$ can be arbitrary. In this paper we limit ourselves to the regime when all three scattering lengths are close to each other, $`a_{11}a_{12}a_{22}`$. In the case when the two species are rubidium atoms in different hyperfine states, these lengths were found experimentally to differ by no more than a few percent. The assumption that $`g_{11}g_{12}g_{22}`$ also introduces considerable technical simplifications in our treatment. We introduce the “average” scattering amplitude

$$g\frac{g_{11}+g_{22}}{2},$$
(10)
and the deviations from the average

$$\delta g(g_{12}g),\delta g^{}g_{11}g=(g_{22}g),$$
(11)
so that $`\delta g\delta g^{}g`$. The stability condition (9) implies that $`\delta g>0`$. Analogously, we introduce the average scattering lengths $`a`$ and the deviations $`\delta a`$ and $`\delta a^{}`$. Note that in the limit $`g_{11}=g_{12}=g_{22}`$ the Hamiltonian has an SU(2) symmetry which leads to interesting implications .

With the Lagrangian at hand, the discussion of symmetry in the Introduction can be made concrete. In the absence of the coupling drive, $`\mathrm{\Omega }=0`$, the Lagrangian (1) possesses a U(1) $`\times `$ U(1) symmetry with respect to independent phase rotations of the fields,

$$\psi _1e^{i\alpha _1}\psi _1,\psi _2e^{i\alpha _2}\psi _2.$$
(12)
The corresponding conservation laws are those of the numbers of particles of each species, $`N_1=d^3𝒙\psi _1^{}\psi _1`$ and $`N_2=d^3𝒙\psi _2^{}\psi _2`$. That $`N_1`$ and $`N_2`$ at $`\mathrm{\Omega }=0`$ are conserved separately is actually a basic assumption made when we wrote down the Lagrangian (1). This assumption is not automatically satisfied: it requires that only elastic scattering between atoms is allowed; inelastic scattering is forbidden. For binary Bose systems made of rubidium atoms, this appears to be a good approximation .

Once the coupling drive is turned on ($`\mathrm{\Omega }0`$), the Lagrangian is invariant only under a subset of the original U(1) $`\times `$ U(1) rotations; namely, those which rotate both $`\psi _1`$ and $`\psi _2`$ by the same angle,

$$\psi _1e^{i\alpha }\psi _1,\psi _2e^{i\alpha }\psi _2.$$
(13)
Therefore one of the U(1) symmetries the system enjoyed at $`\mathrm{\Omega }=0`$ is explicitly violated. Applying the Goldstone theorem, we conclude that, at $`\mathrm{\Omega }=0`$, there are two gapless excitations and only one of these modes remain gapless at $`\mathrm{\Omega }0`$. The gapless modes at $`\mathrm{\Omega }=0`$ are the phonons of the two types of sounds. One corresponds to the ordinary density wave (B mode in our paper, see below), and another to the concentration wave (A mode) in which the densities of the two species oscillate relative to each other in such a way that the total density remains constant.In the hydrodynamic regime these two modes would correspond to the first and the third sounds respectively, according to the established classification . In the zero-temperature limit the density wave is termed Bogolyubov (zero) sound, and there is no established term for the concentration wave. If we view the two components as being made of the same atoms, but in different hyperfine levels, then the mode A can be alternatively interpreted as a spin polarization wave, in which the density of nuclear spin is a function of space and time.

When the coupling drive is on, only the density wave remains gapless; the phonon of the concentration (spin density) wave is gapped.

Let us compute the sound speeds at $`\mathrm{\Omega }=0`$. To this end we write $`\psi _1`$ and $`\psi _2`$ as

$$\psi _1=\sqrt{n_1+\delta n_1}e^{i\phi _1},\psi _2=\sqrt{n_2+\delta n_2}e^{i\phi _2},$$
(14)
and expand Eq. (1) to second order of $`\delta n`$ and $`\phi `$ (we will see that $`\delta n\phi `$). We find

$`L^{(2)}`$ $`=`$ $`\mathrm{}(\delta n_1_t\phi _1+\delta n_2_t\phi _2){\displaystyle \frac{\mathrm{}^2}{2m}}(n_1|\phi _1|^2+n_2|\phi _2|^2)`$ (16)
$`{\displaystyle \frac{1}{2}}(g_{11}\delta n_1^2+2g_{12}\delta n_1\delta n_2+g_{22}\delta n_2^2)`$
(we have thrown away total derivatives). The density fluctuations $`\delta n_i`$ can be “integrated out” and replaced by the saddle point values $`\delta n_i=\mathrm{}_j(g^1)_{ij}_t\phi _j`$. As a result, Eq. (16) becomes

$$L^{(2)}=\frac{\mathrm{}^2}{2}\underset{ij}{}(g^1)_{ij}_t\phi _i_t\phi _j\frac{\mathrm{}^2}{2m}(n_1|\phi _1|^2+n_2|\phi _2|^2).$$
(17)
Thus, the dispersion relations for the phonons are linear, $`\omega =uk`$, and the sound speed $`u`$ satisfies the characteristic equation

$$det\left(\left(g^1\right)_{ij}\frac{1}{mu^2}\left(\begin{array}{cc}n_1& 0\\ 0& n_2\end{array}\right)\right)=0.$$
(18)
When $`g_{11}g_{12}g_{22}`$, the solutions are

$$u_\mathrm{A}^2=\frac{2\delta g}{m}\frac{n_1n_2}{n},u_\mathrm{B}^2=\frac{gn}{m}.$$
(19)
We see that when $`\delta gg`$ the speed of the concentration wave (A) is much smaller than that of the density wave (B), $`u_\mathrm{A}u_\mathrm{B}`$. The system is “stiffer” towards density fluctuations than towards fluctuations of concentration. That modes A and B are indeed concentration and density fluctuations respectively is seen from the corresponding eigenvectors. The A mode corresponds to such fluctuations in which

$$\text{A:}\left(\begin{array}{c}\phi _1\\ \phi _2\end{array}\right)\left(\begin{array}{c}n_2\\ n_1\end{array}\right),\left(\begin{array}{c}\delta n_1\\ \delta n_2\end{array}\right)\left(\begin{array}{c}1\\ 1\end{array}\right),$$
(20)
while the B mode corresponds to

$$\text{B:}\left(\begin{array}{c}\phi _1\\ \phi _2\end{array}\right)\left(\begin{array}{c}1\\ 1\end{array}\right),\left(\begin{array}{c}\delta n_1\\ \delta n_2\end{array}\right)\left(\begin{array}{c}n_1\\ n_2\end{array}\right).$$
(21)
Therefore, in the A sound $`n_1`$ and $`n_2`$ fluctuate in such a way that the overall density remains constant ($`\delta n_1+\delta n_2=0`$), while the B sound corresponds to density waves in which $`n_1/n_2`$, or concentration, is unchanged ($`\delta n_1/n_1=\delta n_2/n_2`$). The Lagrangian (17), in terms of the normal modes

$`\phi _\mathrm{A}\phi _1\phi _2\text{and}\phi _\mathrm{B}(2/n)(n_1\phi _1+n_2\phi _2),`$ (22)
has the form

$$L^{(2)}=\frac{\mathrm{}^2n}{8m}\left\{\frac{4n_1n_2}{n^2}\left[u_\mathrm{A}^2(_t\phi _\mathrm{A})^2(\phi _\mathrm{A})^2\right]+\left[u_\mathrm{B}^2(_t\phi _\mathrm{B})^2(\phi _\mathrm{B})^2\right]\right\}.$$
(23)
When the coupling drive is on, $`\mathrm{\Omega }0`$, one should add to Eq. (17) the potential energy term $`\mathrm{}\mathrm{\Omega }\sqrt{n_1n_2}\mathrm{cos}(\phi _1\phi _2)`$. In the Hamiltonian approach, the effect of the coupling on the spectrum can be treated by diagonalizing Hamiltonian using an exact transformation .  The B sound is not affected by this term since it corresponds to such fluctuations in which $`\phi _1=\phi _2`$. The phonons of the A sound acquire a gap

$$\mathrm{\Delta }=\sqrt{2\mathrm{}\mathrm{\Omega }\delta gn}\left(\frac{n_1n_2}{n^2}\right)^{1/4}$$
(24)
for small values of $`\mathrm{\Omega }`$.

In our further discussion we shall need the formulas for the healing, or correlation, lengths, which are defined via the response of the system to a static source coupled locally to the particle density, $`\mu _1(𝒙)|\psi _1|^2+\mu _2(𝒙)|\psi _2|^2`$. As in the case of the sounds, there are also two healing lengths,

$`\xi _\mathrm{A}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{8}}}{\displaystyle \frac{\mathrm{}}{\sqrt{m\delta gn}}}\left({\displaystyle \frac{n^2}{n_1n_2}}\right)^{1/2}={\displaystyle \frac{1}{\sqrt{2}}}\left({\displaystyle \frac{g}{\delta g}}\right)^{1/2}\left({\displaystyle \frac{n^2}{n_1n_2}}\right)^{1/2}\xi _\mathrm{B},`$ (25)
$`\xi _\mathrm{B}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{\mathrm{}}{\sqrt{mgn}}}={\displaystyle \frac{1}{\sqrt{16\pi an}}}.`$ (26)
As seen from Eq. (26), $`\xi _\mathrm{A}\xi _\mathrm{B}`$. This is because $`\xi _\mathrm{B}`$ is the correlation length of fluctuations of the overall density, while $`\xi _\mathrm{A}`$ is the correlation length of the relative density.

If we take the values typical for present day experiments: $`n10^{14}`$ cm-3, $`a50`$ Å, and $`\delta g/g10^2`$, then $`\xi _\mathrm{B}0.2`$ $`\mu `$m, and $`\xi _\mathrm{A}3`$ $`\mu `$m. These lengths are smaller than the typical system size in experiments.

## III The domain wall of relative phase

The existence of the domain wall of relative phase in binary BEC with a coupling drive can be shown in a rather simple way. Let us first focus only on fluctuations of the fields on length scales much larger than the largest healing length $`\xi _\mathrm{A}`$. In this case, the amplitudes $`n_1`$ and $`n_2`$ of $`\psi _{1,2}`$ can be regarded as “frozen” and the only important degrees of freedom are the phases, i.e., $`\phi _1`$ and $`\phi _2`$. The energy of the system is a functional of $`\phi _1`$ and $`\phi _2`$,

$$E[\phi _1,\phi _2]=d^3𝒙\left[\frac{\mathrm{}^2}{2m}(n_1(\phi _1)^2+n_2(\phi _2)^2)\mathrm{}\mathrm{\Omega }\sqrt{n_1n_2}\mathrm{cos}(\phi _1\phi _2)\right].$$
(27)
The potential energy term $`\mathrm{}\mathrm{\Omega }\sqrt{n_1n_2}\mathrm{cos}(\phi _1\phi _2)`$ has its minimum at $`\phi _\mathrm{A}\phi _1\phi _2=0`$. The configuration in which $`\phi _\mathrm{A}=0`$ over the whole space is obviously the global minimum of the total energy, and hence is the ground state.

The domain wall solution that we shall describe is, in contrast, a local minimum of the energy. To find the profile of the domain wall, we vary Eq. (27) with respect to $`\phi _1`$ and $`\phi _2`$ and obtain the following equations

$$\frac{\mathrm{}^2}{m}n_1^2\phi _1=\frac{\mathrm{}^2}{m}n_2^2\phi _2=\mathrm{}\mathrm{\Omega }\sqrt{n_1n_2}\mathrm{sin}(\phi _1\phi _2).$$
(28)
For domain walls, $`\phi _1`$ and $`\phi _2`$ are functions of only one coordinate, say, $`z`$, and at $`z=\pm \mathrm{}`$ both approach constant values. A nontrivial solution to Eq. (28) satisfying these conditions is

$`\phi _1`$ $`=`$ $`{\displaystyle \frac{n_2}{n}}\phi _\mathrm{A},\phi _2={\displaystyle \frac{n_1}{n}}\phi _\mathrm{A},`$ (29)
$`\text{where}\phi _\mathrm{A}`$ $`=`$ $`4\mathrm{arctan}e^{kz},k^2={\displaystyle \frac{m\mathrm{\Omega }}{\mathrm{}}}{\displaystyle \frac{n}{\sqrt{n_1n_2}}}.`$ (30)
The characteristic width of the domain wall is $`k^1`$. The tension of the domain wall is

$$\sigma =8\frac{\mathrm{}^{3/2}\mathrm{\Omega }^{1/2}n}{m^{1/2}}\left(\frac{n_1n_2}{n^2}\right)^{3/4}.$$
(31)
The relative phase $`\phi _\mathrm{A}`$ changes from 0 to $`2\pi `$ as $`z`$ runs from $`\mathrm{}`$ to $`+\mathrm{}`$. Note that $`\phi _\mathrm{A}`$ is defined modulo $`2\pi `$ so it goes a full circle as one passes through the wall. Therefore from the point of view of the energy functional (27) the domain wall (30) is a topologically nontrivial configuration, which can not be continuously deformed into the $`\phi _\mathrm{A}=0`$ configuration. In fact, one can prove that the domain wall is a configuration with minimal energy defined in Eq. (27) among those where $`\phi _\mathrm{A}`$ changes by $`2\pi `$ from $`z=\mathrm{}`$ to $`z=+\mathrm{}`$, and hence, cannot decay away, as long as Eq. (27) applies. The domain wall is similar to the soliton of the sine-Gordon model. There is a small difference: the ground states on two sides of the domain wall are different: at $`z=\mathrm{}`$ $`\phi _1=\phi _2=0`$, while the ground state at $`z=+\mathrm{}`$ is $`\phi _1=2\pi n_2/n`$, $`\phi _2=2\pi n_1/n`$. Since $`\phi _{1,2}`$ are defined mod $`2\pi `$ the latter is equivalent to $`\phi _1=\phi _2=2\pi n_2/n`$.

In reality, Eq. (27) is not the full Hamiltonian of the system: it is only an effective description valid at length scales larger than both healing lengths. The full theory (3) contains, besides $`\phi _1`$ and $`\phi _2`$, also the density fluctuations $`\delta n_1`$ and $`\delta n_2`$. As a consequence, in the full theory, the domain wall (30) can be continuously deformed into the trivial configuration $`\phi _\mathrm{A}=0`$. Such deformations necessarily pass through field configurations where either $`n_1`$ or $`n_2`$ vanishes at some points: at these points $`\phi _\mathrm{A}`$ is ill-defined. Thus, the domain wall is not truly topological and can “unwind”, i.e., decay away. The fact that the ground states on the two sides of the wall are different does not prevent the decay: it is possible to construct a field configuration interpolating between $`\phi _1=\phi _2=0`$ and $`\phi _1=\phi _2=2\pi n_2/n`$ with arbitrarily small energy per unit area.

Although the domain wall is not a global minimum of the energy functional, it can still be a local minimum. In this case, to deform the domain wall into a “topologically” trivial configuration with $`\phi _\mathrm{A}=0`$ one has to overcome an energy barrier. The wall is in this case metastable. From our previous discussion one can conclude that, roughly speaking, the wall is metastable when Eq. (27) applies and is unstable when Eq. (27) is not applicable. For Eq. (27) to be valid the wall has to be wider than the largest healing length, $`\xi _\mathrm{A}`$. Since the width of the wall decreases as one increases the Rabi frequency $`\mathrm{\Omega }`$, the wall is metastable only when $`\mathrm{\Omega }`$ is smaller than some critical value $`\mathrm{\Omega }_c`$. Let us define $`\mathrm{\Omega }_0`$ as the value of the Rabi frequency at which the width of the wall $`k^1`$, as defined in Eq. (30), is equal to to the longer healing length $`\xi _\mathrm{A}`$ in Eq. (26):

$$\mathrm{}\mathrm{\Omega }_0=8\delta gn\left(\frac{n_1n_2}{n^2}\right)^{3/2}.$$
(32)
The wall is metastable when $`\mathrm{\Omega }`$ is less than some critical value $`\mathrm{\Omega }_c`$ of order $`\mathrm{\Omega }_0`$,

$$\mathrm{\Omega }_c\mathrm{\Omega }_0.$$
(33)
To find the exact value of $`\mathrm{\Omega }_c`$ one needs to perform a more refined stability analysis. We present such an analysis for the special case $`n_1=n_2`$ (i.e., when the densities of the two species are equal) in Appendix A. Parametrically, the result is consistent with Eq. (33), and the ratio $`\mathrm{\Omega }_c/\mathrm{\Omega }_0`$ is found to be $`1/3`$.

Using the numerical values typical for experiments with Rb: $`n`$ a few $`10^{14}`$ cm-3, $`\delta a1`$ Å, one finds $`\mathrm{\Omega }_0100`$ Hz. Therefore to have a stable wall the Rabi frequency needs to be smaller than about 100 Hz. The domain wall cannot be thinner than the longer correlation length $`\xi _\mathrm{A}`$, which was estimated above to be a few $`\mu `$m.

Even when $`\mathrm{\Omega }<\mathrm{\Omega }_c`$, a metastable wall can still spontaneously decay (burst). Such a decay, as we have said, requires overcoming a potential barrier. At sufficiently large temperature, the mechanism for the decay is thermal activation (over the barrier). At zero or small temperatures, the mechanism of the decay is quantum tunneling under the energy barrier. The decay through thermal activation, which is relevant for the temperatures achieved in current experiments, is considered in Appendix B, where it is shown that the decay rate is exponentially suppressed, since a global (macroscopic) fluctuation is required, unless $`\mathrm{\Omega }`$ is very close to $`\mathrm{\Omega }_c`$.

## IV Physical interpretation of the domain wall

The domain wall solution found in Sec. III allows an interesting physical interpretation which suggests a possible way for their creation in experiments . We first note that in a BEC the superfluid velocity is proportional to the gradient of the phase. In a two-component BEC, there are two such velocities.

$$𝒗_1=\frac{\mathrm{}}{m}\mathbf{}\phi _1,𝒗_2=\frac{\mathrm{}}{m}\mathbf{}\phi _2.$$
(34)
Equation (28), in particular, implies that the total particle number current vanishes, $`n_1𝒗_1+n_2𝒗_2=0`$. Individually, however, the particle number current of each species $`𝒋_1n_1𝒗_1`$ and $`𝒋_2n_2𝒗_2`$ are nonzero. From Eq. (30) we find

$$𝒗_1=\frac{n_2}{n}\frac{2k}{\mathrm{cosh}kz}\widehat{𝒛},𝒗_2=\frac{n_1}{n}\frac{2k}{\mathrm{cosh}kz}\widehat{𝒛}.$$
(35)
Thus the domain wall is a configuration where the two components flow in opposite directions. The flow is illustrated in Fig. 1. The velocities of the components are largest at the center of the wall ($`z=0`$) and decrease as one moves toward the edge of the wall. The flow is concentrated on the wall; outside the wall ($`|z|k^1`$) there is essentially no flow.

Equation (28) can be rewritten in terms of the currents as

$$\mathbf{}𝒋_1=\mathbf{}𝒋_2=\mathrm{\Omega }\sqrt{n_1n_2}\mathrm{sin}\phi _\mathrm{A}.$$
(36)
For a stationary configuration as the one we are considering, Eq. (36) means that the number of particles in each species is not conserved. It also implies that there is a conversion of atoms between the two energy levels due to the coupling drive. In the left half of the wall $`\phi _\mathrm{A}<0`$ and there is a conversion of atoms of the second type to atoms of the first type. In the right half $`\phi _\mathrm{A}>0`$ and the conversion goes the opposite way (Fig. 1). The rate of conversion is

$$\mathrm{\Omega }\sqrt{n_1n_2}\mathrm{sin}\phi _\mathrm{A}=2\mathrm{\Omega }\sqrt{n_1n_2}.\frac{\mathrm{sinh}kz}{\mathrm{cosh}^2kz}.$$
(37)
As is the flow, the conversion rate is also maximal near the wall. Far from the wall ($`|z|k^1`$) there is essentially no conversion. The conversion rate (37) changes sign at $`z=0`$. Since different species correspond to different energy levels of the atom, energy is absorbed in one half of the wall and released in the other. §§§ Recalling that, in current experiments, atoms 1 and 2 correspond to different nuclear spin states of the same atom, we can view the wall as a stationary configuration with a flow of (nuclear) spin across $`z=0`$ plane accompanied by spin flips on both sides of it.

The interpretation of the domain walls given above suggests a possible method for their creation in experiments . One starts with the coupling drive off ($`\mathrm{\Omega }=0`$) and prepares a state where the two condensates flow in opposite directions (for example, by manipulating the traps). In such a state the relative phase $`\phi _\mathrm{A}=\phi _1\phi _2`$ is a linear function of the coordinate along the direction of motion (say, $`z`$). One then slowly increases $`\mathrm{\Omega }`$. The domain walls will be created and the centers of the walls are located at the points where $`\phi _\mathrm{A}`$ was an odd multiple of $`\pi `$ ($`\pm \pi `$, $`\pm 3\pi `$, $`\pm 5\pi `$, etc.) before $`\mathrm{\Omega }`$ was turned on. By changing the velocity of the initial relative motion of the condensates and the final values of $`\mathrm{\Omega }`$ one can change the separation between the domain walls and their width. Such controlled creation of the domain walls, hopefully, can be achieved in future experiments.

## V Vortex as the boundary of the wall and vortex confinement

So far we have always considered infinite domain walls which have no boundary. It is also interesting, and perhaps more realistic, to consider domain walls with a boundary. We shall show that the domain wall can be bounded by a vortex line.

Suppose we have a finite domain wall whose boundary is a closed contour $`𝒞`$ (Fig. 2). We shall assume that the length of $`𝒞`$ is much larger than the width of the wall $`k^1`$ so one can view the domain wall as an infinitely thin membrane stretched on $`𝒞`$ (we shall call this picture the “thin-wall approximation”). Let us now take another, smaller contour which has a nontrivial linking with $`𝒞`$ ($`𝒟`$ in Fig. 2). As one goes along $`𝒟`$, one crosses the membrane once, so the relative phase $`\phi _\mathrm{A}`$ changes by $`2\pi `$. This is exactly what one expects from a vortex. Therefore, $`𝒞`$ can be a vortex line. We recall that the size of the core of the vortex is the healing length $`\xi _\mathrm{A}`$, which is smaller or of the same order as the width of the wall. Therefore, in the thin wall approximation, we have an infinitely thin membrane bounded by an infinitely thin vortex line. One should note that such a bounded domain wall will tend to shrink to reduce its energy, which comes from the wall tension and the tension of the boundary vortex.

Conversely, for small non-zero $`\mathrm{\Omega }`$, a vortex must have a domain wall attached to it to minimize the energy due to the nontrivial phase $`\phi _\mathrm{A}`$ winding around it. This means that the energy of a single vortex per unit length is increasing linearly with the size of the system in the transverse direction. This is in contrast to the situation at $`\mathrm{\Omega }=0`$ when the vortex tension has only a logarithmic dependence on the size of the system.

Note that there are two types of vortices in binary BEC. Those of the first type, which we shall call the $`\phi _1`$ vortices, have the condensate $`\psi _1`$ vanishing at the vortex center, while $`\psi _2`$ is nonzero. Analogously the $`\phi _2`$ vortices have $`\psi _2=0`$ and $`\psi _10`$ at their centers. As one goes around a $`\phi _1`$ vortex, $`\phi _1`$ changes by $`2\pi `$, while $`\phi _2`$ does not change, and vice versa for a $`\phi _2`$ vortex. Thus $`\phi _\mathrm{A}`$ changes by either $`2\pi `$ or $`2\pi `$ for the two types of vortices, so the domain wall can be bounded by a vortex of either type.

In contrast to an individual vortex, a pair of $`\phi _1`$ and $`\phi _2`$ vortices, placed parallel to each other, will have energy per unit length which is only logarithmically divergent. That is because the $`\phi _\mathrm{A}`$ “charges” of the two vortices cancel each other, so $`\phi _\mathrm{A}`$ is trivial at spatial infinity (no winding). The same situation occurs for a pair of parallel vortices of the same type ($`\phi _1`$ or $`\phi _2`$) with opposite winding (however, such vortex-antivortex pair can annihilate, while a $`\phi _1\phi _2`$ pair cannot). In a certain sense, one can talk about the phenomenon of “vortex confinement”: vortices exist only in pairs. This confinement should, in principle, be observable experimentally, in a rotating two-component BEC, which can be already created in a laboratory . The vortices are usually identified by the density depletion at their cores, but can be also seen as dislocations of interference fringes due to phase singularities . With the coupling drive off ($`\mathrm{\Omega }=0`$) such a system contains an equal number of $`\phi _1`$ and $`\phi _2`$ vortices which are distributed in space with no particular correlation between $`\phi _1`$ and $`\phi _2`$ vortices. As one turns on $`\mathrm{\Omega }`$, the vortices will start to pair up and at some point the system will become a collection of composite objects, each being a bound state of a $`\phi _1`$ and a $`\phi _2`$ vortex (Fig. 3).

The phenomenon of vortex confinement is very similar to that of quark confinement in the theory of strong interaction (quantum Chromodynamics). Similarly to our vortices, quarks and antiquarks do not exist as individual objects, but are confined into composite objects — hadrons. The analogy with quantum Chromodynamics actually stretches further. If one places a $`\phi _1`$ vortex and a $`\phi _2`$ vortex at a distance larger than $`k^1`$, then a domain wall that connects these two vortices will be formed. The tension of the domain wall is the force, per unit length, that attracts the two vortices. The attractive force between the two vortices is thus independent of their separation, given that the latter is larger than the width of the domain wall. This is analogous to the confining force between a quark and an antiquark, which is also constant at large distances. A confinement model which resembles most the confinement of the vortices is the three-dimensional compact quantum electrodynamics considered in Ref. . In this theory the worldlines of electrically charged particles are analogous to the vortices in BEC.

One can also imagine a system of two vortices which rotate around each other so that the confining force (from the domain wall) balances the centrifugal force. Such a system is analogous to the high-spin meson states in hadronic physics where a quark and an antiquark rotate around each other.

## VI Conclusion

We have shown that in a system of two interpenetrating BEC with a coupling drive, there exists a domain wall solution. The relative phase between the two condensates changes by $`2\pi `$ as one travels through the wall. The wall solution is formally similar to the kink in the sine-Gordon field theory, yet it is not topologically stable and can decay. In this respect, the wall is more similar to a soap film, which can spontaneously burst.

From the mathematical point of view, the domain walls discussed in this paper are similar to the ones which have been studied in particle physics. Such domain walls appear at least in two contexts: in the theory of the hypothetical axion  and in high-density quark matter .Another case discussed in particle-physics literature, domain walls in zero-density QCD , is beyond theoretical control. The similarity is that the domain wall solution arises from the spontaneous breaking of an approximate U(1) symmetry. In all cases the domain wall exists only if the explicit violation of the U(1) symmetry (determined by the value of the Rabi frequency in our case) is small enough. The decay of the wall in all the examples occurs via hole nucleation.

As to the experimental realization of the domain wall, one could be optimistic since the estimated critical value of the Rabi frequency, for densities and scattering lengths typical for the rubidium gas in recent experiments, is of the order 100 Hz which is not too small. The width of the wall, which might be as small as a few $`\mu `$m if the Rabi frequency is not much smaller than critical, can be also accommodated inside condensates of the size characteristic of present-day experiments. The apparent immediate obstacle is still the creation of a system where two BEC truly interpenetrate. In the recent experiments with trapped atomic gas the region of overlap between the two BEC is still small. One can hope the technical problems of making a genuine two-component BEC to be solved in the near future which would enable one to study the domain walls experimentally.

###### Acknowledgements.

The authors thank Richard Friedberg for discussions. We thank RIKEN, Brookhaven National Laboratory, and U.S. Department of Energy \[DE-AC02-98CH10886\] for providing the facilities essential for the completion of this work. The work of DTS is supported, in part, by a DOE OJI grant.

## A Wall’s local stability, $`n_1=n_2`$

In the case when the densities of the two components are equal, $`n_1=n_2`$, the maximal frequency $`\mathrm{\Omega }_c`$, below which the wall configuration is still locally stable, can be found analytically. For illustrative purposes we shall consider this particular case in details.

Let us recall that when the scattering lengths $`a_{11}`$, $`a_{12}`$, and $`a_{22}`$ are approximately equal (assuming $`a_{11}a_{22}a_{12}^2>0`$) there are two healing lengths $`\xi _\mathrm{A}`$ and $`\xi _\mathrm{B}`$. The healing length related to fluctuations of the overall density $`\xi _\mathrm{B}`$ is much smaller than the one related to fluctuations of the relative density $`\xi _\mathrm{A}`$: $`\xi _\mathrm{B}\xi _\mathrm{A}`$. We shall be interested in the case when the width of the wall $`k^1`$ is much larger than $`\xi _\mathrm{B}`$ (but we shall not presume any relation between $`k^1`$ and $`\xi _\mathrm{A}`$). In this case, the total density $`n=|\psi _1|^2+|\psi _2|^2`$ can be considered as frozen, and the system can be described in terms of three variables: $`\theta `$, $`\phi _1`$ and $`\phi _2`$,

$$\left(\begin{array}{c}\psi _1\\ \psi _2\end{array}\right)=\sqrt{n}\left(\begin{array}{c}\mathrm{cos}(\theta /2)e^{i\phi _1}\\ \mathrm{sin}(\theta /2)e^{i\phi _2}\end{array}\right),$$
(A1)
where $`\theta `$ runs from $`0`$ and $`\pi `$. In terms of these three variables, the Lagrangian has the form

$$L=\mathrm{}n[\mathrm{cos}^2(\theta /2)_t\phi _1+\mathrm{sin}^2(\theta /2)_t\phi _2]H,$$
(A2)
where $`H`$ is the following Hamiltonian

$`H(\theta ,\phi _1,\phi _2)`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2n}{2m}}\left[{\displaystyle \frac{1}{4}}(\theta )^2+\mathrm{cos}^2{\displaystyle \frac{\theta }{2}}(\phi _1)^2+\mathrm{sin}^2{\displaystyle \frac{\theta }{2}}(\phi _2)^2\right]+{\displaystyle \frac{1}{2}}(\delta \mu n\delta g^{}n^2)\mathrm{cos}\theta `$ (A4)
$`{\displaystyle \frac{1}{4}}\delta gn^2\mathrm{sin}^2\theta {\displaystyle \frac{1}{2}}\mathrm{}\mathrm{\Omega }n\mathrm{sin}\theta \mathrm{cos}(\phi _1\phi _2).`$
In Eq. (A4) $`\delta \mu =\mu _1\mu _2`$. The ground state is found by minimizing the potential term in Eq. (A4). From here on we consider the special case when in the ground state the density of atoms of the two species are equal, $`n_1=n_2=n/2`$ (or equivalently $`\theta =\pi /2`$.) This requires $`\delta \mu =\delta g^{}n`$. In this case $`H`$ possesses a discrete symmetry with respect to replacing $`\theta \pi \theta `$, $`\phi _1\phi _2`$. This symmetry is what makes it possible to find $`\mathrm{\Omega }_c`$ analytically.

It is more convenient to use, instead of $`\phi _{1,2}`$, the normal modes $`\phi _{\mathrm{A},\mathrm{B}}`$ defined in Eqs. (22), which in the case $`n_1=n_2=n/2`$ have the form

$`\phi _\mathrm{A}`$ $`=`$ $`\phi _1\phi _2,`$ (A5)
$`\phi _\mathrm{B}`$ $`=`$ $`\phi _1+\phi _2.`$ (A6)
In terms of these variables the Hamiltonian becomes

$`H`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2n}{8m}}\left[(\theta )^2+(\phi _\mathrm{A})^2+(\phi _\mathrm{B})^2+2\mathrm{cos}\theta \mathbf{}\phi _\mathrm{A}\mathbf{}\phi _\mathrm{B}\right]`$ (A8)
$`{\displaystyle \frac{1}{4}}\delta gn^2\mathrm{sin}^2\theta {\displaystyle \frac{1}{2}}\mathrm{}\mathrm{\Omega }n\mathrm{sin}\theta \mathrm{cos}\phi _\mathrm{A}.`$
In order to find the domain wall configuration, we need to extremize the energy with respect to $`\theta `$ and $`\phi _{\mathrm{A},\mathrm{B}}`$. Varying with respect to $`\phi _\mathrm{B}`$, one finds

$$\mathbf{}(\mathbf{}\phi _\mathrm{B}+\mathrm{cos}\theta \mathbf{}\phi _\mathrm{A})=0.$$
(A9)
This equation determines $`\phi _\mathrm{B}`$ for given $`\theta `$ and $`\phi _\mathrm{A}`$. The task of solving Eq. (A9) becomes much simpler if one assumes that all variables depend only one coordinate $`z`$. In this case

$$_z\phi _\mathrm{B}+\mathrm{cos}\theta _z\phi _\mathrm{A}=0,$$
(A10)
which can be trivially solved: $`\phi _\mathrm{B}=^z𝑑z\mathrm{cos}\theta _z\phi _\mathrm{A}`$. After eliminating $`\phi _\mathrm{B}`$, the energy functional one has to minimize is

$$H=\frac{\mathrm{}^2n}{8m}\left[(\theta )^2+\mathrm{sin}^2\theta (\phi _\mathrm{A})^2\right]\frac{1}{4}\delta gn^2\mathrm{sin}^2\theta \frac{1}{2}\mathrm{}\mathrm{\Omega }n\mathrm{sin}\theta \mathrm{cos}\phi _\mathrm{A}.$$
(A11)
It is easy to check that the following configuration is always a local extremum of Eq. (A11):

$`\overline{\theta }(z)`$ $`=`$ $`{\displaystyle \frac{\pi }{2}},`$ (A12)
$`\overline{\phi }_\mathrm{A}(z)`$ $`=`$ $`4\mathrm{arctan}e^{kz},k^2={\displaystyle \frac{2m\mathrm{\Omega }}{\mathrm{}}}.`$ (A13)
Equation (A12) can be guessed from the symmetry of Eq. (A11) under $`\theta \pi \theta `$. To see if the domain wall solution is a local minimum, one needs to expand $`H`$ in the vicinity of (A12,A13). One writes

$`\theta (z)`$ $`=`$ $`\overline{\theta }(z)+\stackrel{~}{\theta }(z),`$ (A14)
$`\phi _\mathrm{A}(z)`$ $`=`$ $`\overline{\phi }_\mathrm{A}(z)+\stackrel{~}{\phi }_\mathrm{A}(z).`$ (A15)
To the second order in $`\stackrel{~}{\theta }`$ and $`\stackrel{~}{\phi }_\mathrm{A}`$ the Hamiltonian (A11) is

$$H^{(2)}=\frac{\mathrm{}^2n}{8m}\left[(\stackrel{~}{\theta })^2+\left(\frac{2\delta gnm}{\mathrm{}^2}+k^2\mathrm{cos}\overline{\phi }_\mathrm{A}(\overline{\phi }_\mathrm{A})^2\right)\stackrel{~}{\theta }^2\right]+\frac{\mathrm{}^2n}{8m}\left[(\stackrel{~}{\phi }_\mathrm{A})^2+k^2\mathrm{cos}\overline{\phi }_\mathrm{A}\stackrel{~}{\phi }_\mathrm{A}^2\right].$$
(A16)
One has to find the eigenmodes of Eq. (A16): if there are no negative modes then the domain wall is a local minimum of the energy; if there exist a negative mode then the domain wall is unstable. The second, $`\theta `$-independent, term in Eq. (A16) does not have negative modes (it has only one zero mode corresponding to the translation of the wall along the $`z`$ direction) and does not lead to instability, and hence can be ignored. Taking into account the explicit solution $`\overline{\phi }_\mathrm{A}=4\mathrm{arctan}e^{kz}`$, the first term in Eq. (A16) is

$$\frac{\mathrm{}^2n}{8m}\left[(\stackrel{~}{\theta })^2+\left(\frac{2\delta gnm}{\mathrm{}^2}+k^2\frac{6k^2}{\mathrm{cosh}^2kz}\right)\stackrel{~}{\theta }^2\right].$$
(A17)
The well-known operator

$$^2\frac{6k^2}{\mathrm{cosh}^2kz}$$
(A18)
has the lowest eigenvalue equal to $`4k^2`$, corresponding to the eigenfunction $`\mathrm{cosh}^2kz`$, which implies that $`H`$ does not have a negative mode if $`\mathrm{\Omega }<\mathrm{\Omega }_c`$, where

$$\mathrm{\Omega }_c=\frac{1}{3}\frac{\delta gn}{\mathrm{}}.$$
(A19)
When $`\mathrm{\Omega }>\mathrm{\Omega }_c`$, the configuration (A13) is not a local minimum of the energy functional: the domain wall does not exist. The value (A19) is of the same order as $`\mathrm{\Omega }_0`$ in Eq. (32): $`\mathrm{\Omega }_c=\frac{1}{3}\mathrm{\Omega }_0`$. The numerical value $`\frac{1}{3}`$ is specific for the case $`n_1=n_2`$; if $`n_1n_2`$ then $`\mathrm{\Omega }_c/\mathrm{\Omega }_0`$ is, in general, different.

## B Decay of the domain wall

As we have seen above, as long as $`\mathrm{\Omega }<\mathrm{\Omega }_c`$, the wall minimizes the energy of the system with respect to small local variations of the condensates. However, the global minimum of the energy is achieved when the phases of the condensates are constant in space. Since, as discussed above, the wall configuration can be continuously deformed into the global minimum energy configuration, i.e., the wall can decay. Here we shall estimate the lifetime of the wall due to such a decay.

Since the wall minimizes the energy locally, such a deformation necessarily goes through a potential barrier. Thus the decay can only occur by a quantum tunneling through this barrier or, at finite temperature, by a thermal fluctuation over the barrier. In both cases, the decay rate is exponentially suppressed: in the first case by the WKB factor $`e^{S/\mathrm{}}`$, and in the second case by the Boltzmann factor $`e^{E_c/k_BT}`$, where $`E_c`$ is the height of the barrier. The first formula applies at sufficiently low temperature, while the second one applies at higher temperatures. A crude estimate (see below) suggests that the corresponding crossover temperature is quite small (nanokelvins for parameters typical for present-day experiments), and is lower than temperatures achieved in present-day experiments. Thus we shall limit our discussion to the decay by thermal fluctuation. This case is also simpler theoretically. The quantum mechanical decay, relevant for very low temperatures, will be left beyond the scope of this paper.

Assuming the temperature is much smaller than the critical temperature at which one of the condensates melts (so that most atoms are still in the condensates), the rate of thermal activation across the barrier is $`\mathrm{\Gamma }e^{E_c/k_BT}`$, where $`E_c`$ is the height of the barrier at zero temperature. Since we are dealing with an infinite-dimensional configuration space, $`E_c`$ should be understood as the energy at the lowest point of the barrier. This point is a saddle point (the energy has a single negative curvature direction).

Some information about the exponent $`E_c/k_BT`$ can be obtained by a simple scaling argument without a detailed calculation. In terms of the dimensionless variables $`\stackrel{~}{𝒙}`$ defined as

$$𝒙=\frac{\mathrm{}}{\sqrt{m\delta gn}}\stackrel{~}{𝒙},$$
(B1)
the energy becomes

$`E[\theta ,\phi _1,\phi _2]={\displaystyle \frac{\mathrm{}^3n^{1/2}}{m^{3/2}\delta g^{1/2}}}{\displaystyle d^3\stackrel{~}{𝒙}}`$ $`\{{\displaystyle \frac{1}{2}}[{\displaystyle \frac{1}{4}}(\stackrel{~}{}\theta )^2+\mathrm{cos}^2{\displaystyle \frac{\theta }{2}}(\stackrel{~}{}\phi _1)^2+\mathrm{sin}^2{\displaystyle \frac{\theta }{2}}(\stackrel{~}{}\phi _2)^2]`$ (B3)
$`{\displaystyle \frac{1}{4}}(\mathrm{sin}^2\theta +2\mathrm{cos}\theta _0\mathrm{cos}\theta ){\displaystyle \frac{1}{2}}{\displaystyle \frac{\mathrm{}\mathrm{\Omega }}{\delta gn}}\mathrm{sin}\theta \mathrm{cos}(\phi _1\phi _2)\},`$
where $`\theta _0`$ determines the relative condensate densities in the ground state at $`\mathrm{\Omega }=0`$: $`\mathrm{cos}\theta _0=(n_1n_2)/n`$. Consider the dependence on $`\mathrm{\Omega }`$ at a given $`\theta _0`$. The energy functional Eq. (B3) is equal to a dimensionful constant times a dimensionless functional which depends on $`\mathrm{\Omega }`$ only via the ratio $`\mathrm{\Omega }/\mathrm{\Omega }_0`$, where $`\mathrm{\Omega }_0=(\delta gn/\mathrm{})\mathrm{sin}^{3/2}\theta _0`$, same as in (32). Thus the saddle point of the energy is also a function of this dimensionless ratio:

$$E_c=\frac{\mathrm{}^3n^{1/2}}{m^{3/2}\delta g^{1/2}}F\left(\frac{\mathrm{\Omega }}{\mathrm{\Omega }_0}\right).$$
(B4)
If, moreover, we define a temperature $`T_0`$ as

$$k_BT_0=\frac{\mathrm{}^2n^{2/3}}{m},$$
(B5)
which is of the same order as the critical temperature of the Bose-Einstein condensation, then the decay rate for a given $`n_1/n_2`$ (or $`\theta _0`$) can be written as

$$\mathrm{\Gamma }e^{E_c/k_BT}=\mathrm{exp}\left[\frac{\mathrm{}}{m^{1/2}\delta g^{1/2}n^{1/6}}\frac{T_0}{T}F\left(\frac{\mathrm{\Omega }}{\mathrm{\Omega }_0}\right)\right]=\mathrm{exp}\left[\frac{1}{\sqrt{4\pi \delta an^{1/3}}}\frac{T_0}{T}F\left(\frac{\mathrm{\Omega }}{\mathrm{\Omega }_0}\right)\right].$$
(B6)
The form of the function $`F(\mathrm{\Omega }/\mathrm{\Omega }_0)`$ cannot be found from scaling arguments alone. However, when $`\mathrm{\Omega }\mathrm{\Omega }_0`$ one can expect $`F(\mathrm{\Omega }/\mathrm{\Omega }_0)1`$, and then the exponent $`E_c/k_BT`$ is large, since $`\delta an^{1/3}an^{1/3}1`$, and $`TT_0`$.

The function $`F(\mathrm{\Omega }/\mathrm{\Omega }_0)`$ can be computed in the regime $`\mathrm{\Omega }\mathrm{\Omega }_0`$, where the saddle point configuration can be found in the “thin-wall” approximation. In this approximation field configurations are described at length scales much larger than the width of the domain wall. From this point of view the domain wall is an infinitely thin membrane. The saddle point configuration is a membrane with a round hole in it (Fig. 4). The radius of the hole $`R`$ must be much larger than the width of the wall for the thin wall approximation to be valid. As discussed in Sec. V, the rim of the hole must be a vortex, since it is the boundary of the domain wall.

There are two contributions to the energy difference $`E_c`$ between the saddle point and the domain wall configurations. One contribution is negative and comes from the hole (since the hole is the absence of the wall). Another contribution is positive and comes from the rim. Therefore,

$$E=2\pi R\nu \pi R^2\sigma ,$$
(B7)
where $`\nu `$ is the energy per unit length of the vortex (the vortex tension) and $`\sigma `$ is the domain wall tension. The energy (B7) has a maximum when the radius of the hole is

$$R_c\nu /\sigma .$$
(B8)
This is the radius of the critical hole. Indeed, if a hole with a radius $`R>R_c`$ is nucleated, then it will expand and eventually eat up the whole wall. If, in contrast, the radius of the nucleated hole is less than critical, then the hole will shrink and disappear.

Substituting $`R=R_c`$ in Eq. (B7), we find the height of the energy barrier

$$E_c=\frac{\pi \nu ^2}{\sigma }.$$
(B9)
As there are two types of vortices, in Eq. (B9) $`\nu `$ refers to the vortex with the smaller tension. The tension of a straight vortex is logarithmically divergent:

$$\nu _i=\frac{\mathrm{}^2n_i}{2m}2\pi \mathrm{ln}\frac{R}{\xi _\mathrm{A}}.$$
(B10)
where the index $`i=1,2`$ refers to the two types of vortices, and $`R`$ is the long-distance cutoff. The role of $`R`$ is played by either the size of the critical hole $`R_c`$ or the width of the wall $`k^1`$. We shall see at the end of this Section that the two lengths differ only by a logarithm, which does not affect Eq. (B10). Therefore, to logarithmic accuracy, the vortex tension is

$$\nu _i=\frac{\pi \mathrm{}^2n_i}{m}\mathrm{ln}\frac{1}{k\xi _\mathrm{A}}.$$
(B11)
The argument of the logarithm is $`(k\xi _\mathrm{A})^1=(\mathrm{\Omega }_0/\mathrm{\Omega })^{1/2}`$, and is large when $`\mathrm{\Omega }\mathrm{\Omega }_0`$, which justifies the use of the logarithmic approximation.

Without loss of generality, we can assume that $`n_1n_2`$. Then the vortex of the first type has the smallest $`\nu `$ (B10). Substituting Eqs. (B11) and (31) into Eq. (B9), one finds that the barrier height has the form of Eq. (B4), where

$$F\left(\frac{\mathrm{\Omega }}{\mathrm{\Omega }_0}\right)=\frac{\pi ^3}{64\sqrt{2}}\frac{nn_1^{1/2}}{n_2^{3/2}}\left(\frac{\mathrm{\Omega }_0}{\mathrm{\Omega }}\right)^{1/2}\mathrm{ln}^2\frac{\mathrm{\Omega }_0}{\mathrm{\Omega }}.$$
(B12)
The decay rate (B6) can be rewritten in the following form

$$\mathrm{\Gamma }_{\mathrm{thermal}}\mathrm{exp}\left\{\frac{\pi ^2}{128\sqrt{2}}\left[\zeta \left(\frac{3}{2}\right)\right]^{2/3}\left(\frac{n}{n_2}\right)^{3/2}\frac{1}{\sqrt{4\pi \delta an_1^{1/3}}}\frac{T_{c1}}{T}\left(\frac{\mathrm{\Omega }_0}{\mathrm{\Omega }}\right)^{1/2}\mathrm{ln}^2\frac{\mathrm{\Omega }_0}{\mathrm{\Omega }}\right\},$$
(B13)
where we introduced

$$k_BT_{c1}=\frac{2\pi \mathrm{}^2}{m}\left(\frac{n_1}{\zeta (3/2)}\right)^{2/3}$$
(B14)
— the critical temperature for the smallest of the two condensates ($`n_1`$ by our choice). The rate (B13) is exponentially suppressed when $`\mathrm{\Omega }\mathrm{\Omega }_0`$ and $`TT_{c1}`$.

To check the consistency of our assumptions, we note that the radius of the critical hole is

$$R_c=\frac{\nu _1}{\sigma }=\frac{\pi }{16}\sqrt{\frac{\mathrm{}}{m\mathrm{\Omega }}}\frac{n}{n_2}\left(\frac{n_1n_2}{n^2}\right)^{1/4}\mathrm{ln}\frac{\mathrm{\Omega }_0}{\mathrm{\Omega }}.$$
(B15)
Comparing with the width of the wall in Eq. (30), we find that $`R_c`$ is larger than $`k^1`$ by a factor $`\mathrm{ln}(\mathrm{\Omega }_0/\mathrm{\Omega })`$, which is assumed to be parametrically large. Thus, the use of the thin wall approximation is justified.

To make a good estimate of the crossover temperature, below which crossover proceeds via quantum tunneling, rather than via thermal activation, one needs to estimate the action on the tunneling trajectory and compare it to the exponent in the thermal activation rate $`E_c/k_BT`$. For a very crude estimate we can take $`SE_cR_c/u_A`$. The exponents $`S/\mathrm{}`$ and $`E_c/k_BT`$ become comparable at a temperature of order few nanokelvin, given typical parameters of present day experiments.


format: text

text:
# Equilibrium valleys in spin glasses at low temperature

## I Introduction

In the past two decades there has been a large amount of theoretical, experimental, and numerical work on spin glasses (see for example , and references therein). Two different points of view have dominated the frameworks used by most researchers. On one side, there is the idea that mean field features , that are astonishingly different from those encountered in systems without frustration, are relevant in realistic, finite dimensional spin glasses. On the other side, there is the expectation that the (non-trivial) extension  of scaling ideas to these frustrated systems should be the correct framework for short range spin glasses. The results contained in this note will give support to the point of view that many features of the finite dimensional Edwards-Anderson (EA)  spin glasses are common to the ones of the mean field Parisi solution and are well described in terms of the Replica Symmetry Breaking (RSB) approach. This point of view is already supported by a large number of numerical simulations and by studies of ground states. Indeed, numerical simulations indicate, among other features, that the EA model has an average overlap probability distribution $`P(q)`$ which is broad in the thermodynamic limit, signaling (continuous) RSB (see  and references therein). Similarly, ground state computations  show that the energy landscape of the EA model is characterized by system-size valleys whose excitation energies remain of order one in the infinite volume limit. The consequence of this property on the system at temperature $`T>0`$ is likely to be RSB. This is typically expected to imply that for a given instance of the disorder (one given realization of the quenched random couplings) spin configurations will cluster into multiple “valleys”, not related by the up-down symmetry (there is no magnetic field so we shall always work modulo the $`\pm 1`$ symmetry).

The purpose of this work is to study the $`3`$-d EA model at low temperature, and to characterize its valleys, both thermodynamically and for their position-space (i.e., spatial) properties. Among other things, we ask whether valleys having $`O(1)`$ free-energy differences have similar energies and entropies; if not, RSB implies that the fluctuations in these quantities must nearly cancel. Also, what are the typical differences among valleys? How does one go from one valley to another? Before addressing these questions, we have to be far more specific in defining a valley. Our definition is algorithmic, based on a method for clustering configurations according to their spin overlaps (or equivalently according to their Hamming distances). Such a clustering scheme has been used successfully by Hed et al.  for analyzing spin glass ground state configurations. By applying such a clustering algorithm, we determine at a given temperature the “valleys” that characterize each disorder realization. Then we compute various statistical properties of these valleys and of their associated TAP-like states. Interestingly, the picture obtained from looking at $`O(1)`$ energy system-size excitations above the ground state  also applies to our finite temperature valleys. One goes from one valley to another by flipping system-size clusters of spins that are space filling and topologically highly non-trivial: this justifies calling them spongy clusters .

The outline of this paper is as follows. First we specify how our equilibrium configurations were generated. Then, given two such configurations, we consider various ways to compare them with a special attention on topological properties of the set of spins where they differ. We construct the clusters of these spins and classify them into sponge-like, droplet-like, or neither of these. Our data indicate that, if the statistics of such topological events are monotonic with $`L`$, sponge-like events arise with a strictly positive probability in the large $`L`$ limit. Then we describe how we define the valleys into which we will cluster our equilibrium configurations (we give additional details about the clustering algorithm in the appendix). Not surprisingly, configurations assigned to significantly different valleys nearly always have sponge-like differences, i.e., they do not differ in localized, droplet like areas, but in a diffuse set of sites. We compute the (extensive) internal energy of each valley: we find that this quantity fluctuates dramatically from valley to valley, growing very clearly with the lattice size $`L`$. Finally we consider the properties of the TAP-like states associated with the valleys. We find that again the difference among different TAP-like states is typically given by a spongy cluster, just as in the zero temperature case. We also study the link overlap distribution among these states.

## II The model and its numerical simulation

We consider the $`3D`$ Edwards-Anderson (EA) spin glass model, whose Hamiltonian is

$$H=\underset{ij}{}J_{ij}S_iS_j,$$
(1)
where $`ij`$ indicates that the sum is over nearest neighbors. We work on simple cubic lattices with periodic boundary conditions and zero magnetic field. The quenched couplings $`J_{ij}`$ are independent random variables that can take the two values $`\pm 1`$ with equal probability $`1/2`$. In this study we have used spin configurations at $`T=0.5`$. (Note that $`T_c1.1`$.) We will present data taken on lattices of linear size $`L=6`$, $`8`$ and $`12`$.

We rely on a parallel tempering updating scheme (see for example  and references therein) that is very effective for simulating EA spin glasses. The different temperatures used in the runs are uniformly spaced. We use $`19`$ temperature values going from $`T_{min}=0.5`$ to $`T_{MAX}=2.3`$ on the $`L=6`$ lattices, $`49`$ temperature values going from $`T_{min}=0.5`$ to $`T_{MAX}=2.1`$ on the $`L=8`$ lattices, and $`64`$ temperature values going from $`T_{min}=0.5`$ to $`T_{MAX}=2.075`$ on the $`L=12`$ lattices. We have checked in detail that we reach a very good thermalization according to all the usual tests : the acceptance ratio of the tempering sweeps is high (always of order $`0.7`$), the permanence histograms are very flat, and observables look very constant when plotted as a function of the log of the elapsed Monte Carlo time. In addition, we have checked that on individual disorder instances the $`\pm 1`$ symmetry of $`P_J(q)`$ is very well satisfied, see for example figures 1 and 2. In short, since we have used very long MC runs with very conservative and safe choices of the parameters, thermalization seems to be complete.

For each disorder instance and lattice size, we perform the parallel tempering independently on two (sets of) replicas. For each, we first run $`10^6`$ MC sweeps just to reach thermalization; a sweep consists of $`L^3`$ trial spin flips and one tempering $`\beta `$ trial update. Following that, we run for $`1.110^6`$ MC sweeps to gather (equilibrium) statistics. During these sweeps we save $`1100`$ spin configurations, that is one every $`1000`$ MC sweeps, giving a total of $`2200`$ equilibrium configurations for that disorder sample. These are the configurations that are used for the analysis described in the rest of this paper. The runs to generate these configurations took on the order of $`2`$ weeks of machine time on a dedicated Linux cluster based on $`16`$ Pentium II $`450MHz`$ processors.

For each lattice size ($`L=6`$, $`8`$, and $`12`$), we thus have $`16\times 32=512`$ different disorder samples (the factor $`32`$ comes from the use of multi-spin coding on each of the $`16`$ machines), and a total of $`2\times 1100`$ equilibrium spin configurations for each disorder sample.

## III Droplet and sponge differences of equilibrium configurations

Each spin configuration $`C`$ appears at equilibrium, at temperature $`T`$, with probability $`P_B(C)=Z_J^1\mathrm{exp}(\frac{H(C)}{T})`$ (where $`Z_J`$ is the partition function which normalizes the Boltzmann weight for that instance). Consider a pair of spin configurations $`C^{(1)}`$ and $`C^{(2)}`$, selected independently according to their à priori probability, and let $`G`$ be the set of sites where the spins of the two configurations differ. One defines the spin overlap by

$$q(C^{(1)},C^{(2)})=\frac{1}{N}\underset{i=1}{\overset{N}{}}S_i^{(1)}S_i^{(2)},$$
(2)
where in our $`3`$-d case $`N=L^3`$. If $`|G|`$ is the cardinality of $`G`$, then $`q=12|G|/N`$. In figure 1 we show an example (for $`L=12`$) of a distribution $`P_J(q)`$ of overlaps for one realization of the disorder (the dotted curve), and with the solid curve we show the distribution averaged over the disorder, that is over our $`512`$ disorder samples. Because of the up-down symmetry, in the rest of this work we will only look at the $`q0`$ part of the support, and we will ignore the part for negative $`q`$ values. For the disorder instance used in figure 1, most of the overlaps are close to the two values associated with the positions of the two peaks of $`P_J(q)`$. This suggests that in this instance the spin configurations fall into two “valleys”, as will be confirmed later. In figure 2 we consider a different disorder instance where $`P_J(q)`$ has three peaks (always in the $`q>0`$ sector). Because of this kind of behavior, we should investigate what happens when we cluster the configurations according to their mutual overlaps; however, before doing so, let us first discuss how two equilibrium configurations differ spatially.

When two configurations are in the same valley, we expect them to differ only by droplet-like objects; the set of spins $`G`$ forming their difference should consist of small clusters. On the contrary, when considering configurations in two different valleys, we expect $`G`$ to contain a large and topologically non-trivial cluster. To quantify this, we follow reference  and declare $`G`$ to be sponge-like if it and its complement wind around the lattice in all three directions, whereas we call it droplet-like if it does not wind around any of the three directions of the lattice. Finally we label intermediate the other cases of windings. Our data show very clearly a strong correlation between the (relative) size of $`G`$ and its topological class: the larger $`|G|/L^3`$, the more likely it is to be sponge-like. Do sponge-like differences between equilibrium configurations survive in the large $`L`$ limit? In table I we give the probabilities of the three topological classes as a function of the lattice size. Just as in the $`T=0`$ energy landscape studies , the frequencies of sponge-like events increases as $`L`$ increases, suggesting that in the limit $`L\mathrm{}`$ equilibrium configurations differ by sponge-like clusters with a strictly positive probability.

Coming back to the distribution of overlaps, let us decompose $`P(q)`$ into the sum of three densities, one for each topological class:

$$P(q)=\rho _s(q)+\rho _d(q)+\rho _i(q),$$
(3)
where $`\rho _s`$ is the sponge density, $`\rho _d`$ the droplet density, and $`\rho _i`$ the intermediate density. (A different decomposition based on valleys has been introduced by Hed et al. ; our method has the great advantage that it does not rely on any clustering into valleys and is thus parameter-free.) Now consider only those pairs of configurations where $`G`$ is sponge-like: we show in figure 3 the density distribution $`\rho _s`$ of the overlaps belonging to this topological class. The main trend is the widening of the distribution’s support: sponge-like $`G`$’s can occur for increasing values of $`q`$ as $`L`$ grows. Other than that, for small $`q`$, $`\rho _s`$ is stable when $`L`$ increases from $`6`$ to $`12`$: this constitutes further numerical evidence for replica symmetry breaking in the three-dimensional EA model.

## IV Clustering equilibrium configurations

We want to give an appropriate definition of a valley. One would like a valley to be a connected region in phase space that contributes a finite measure to the partition function (this statement can be made precise by coarse graining the phase space). Configurations belonging to the same valley are then expected to be close (i.e., they have a large overlap), while configurations belonging to different valleys are expected to be far (small overlap). Now, in more precise terms, given a collection of equilibrium configurations, we seek to cluster them into families which hereafter we shall refer to as valleys. Although this kind of classification problem is generally ill-posed, one expects good clustering methods to lead to nearly identical results when the valleys are well separated. Such an “ideal” case arises in models with one-step RSB: there, the overlap between two equilibrium configurations converges in the thermodynamic limit to one of two possible values, so assigning configurations to clusters is straight-forward.

In the case of continuous RSB, the situation is slightly ambiguous because in the thermodynamic limit we expect that there will be some hierarchical organization of the states: there will be valleys inside valleys inside valleys… Nevertheless, it is still meaningful to define valleys, as long as the splitting of valleys into sub-valleys is not relevant for the observables considered. If that is indeed the case, one can in fact allow for valley sub-divisions ad infinitum as was done by Hed et al. : those authors clustered ground state configurations in the $`3`$-d $`\pm J`$ EA model, letting the sub-valley sizes become arbitrarily small. Since our goal here is to focus on well separated valleys, we construct our valleys up to some cut-off overlap value $`q^{}`$; in effect, all sub-valleys with overlaps larger than $`q^{}`$ are lumped together. After, we will check by using different values of $`q^{}`$ that our findings are not sensitive to that cut-off; this can be so of course only when considering observables associated with overlaps smaller than that cut-off.

For completeness, we give here a high-level description of our clustering method. The appendix gives some details about the algorithm itself (a similar approach was defined by Iba and Hukushima ). Conceptually, our assignment of equilibrium configurations to a given valley is motivated by the TAP states arising in the mean field picture. Each of our valleys will be associated with a “TAP-like” state which for our purposes is just a set of magnetizations, one for each site. In analogy with TAP states, these magnetizations are the thermal averages of the spin variables, restricted to configurations belonging to the valley of interest; this is like using a “finite volume equilibrium state” to define the TAP states. Thus, given the set of configurations assigned to a valley, we compute the magnetization $`m_i`$ at site $`i`$ by finding the mean of the values of the $`S_i`$ over these configurations. If we want to have $`M`$ valleys, we are to find $`M`$ TAP-like states; algorithmically, we find these $`M`$ states and their associated valleys by iteration. The self-consistency conditions are that:

1. the TAP-like state of valley $`k`$ is defined by the magnetizations $`m_i`$ of that valley;
2. a configuration belongs to the valley $`k`$ if the TAP-like state nearest to it is number $`k`$ (the nearness is measured by the overlap) and if that overlap is greater than $`q^{}`$.

If $`q^{}`$ is close to $`1`$ only very similar configurations will be included into the same valley; on the contrary, if $`q^{}`$ is small most of the time all the configurations will be assigned to just one valley. Intuitively, one can think of $`q^{}`$ as being a resolution; as one increases $`q^{}`$, valleys get resolved into sub-valleys, and as $`q^{}1`$, one reaches the limit where each valley corresponds to a single spin configuration.

Our self-consistent procedure does not generally determine the clustering uniquely, but as was said before, any reasonable clustering algorithm should do well if the valleys are sufficiently separated (that is the case we want to consider in this work). In practice we find the valleys (and their TAP-like states) iteratively: we start by defining a first valley, then when possible and needed we add and define a second valley, and so on. For a given value of the cut-off $`q^{}`$, one clusterizes more and more of the configurations when increasing the number of TAP-like states. One stops the construction when any of the following conditions are met: (i) all spin configurations have been assigned to valleys; (ii) a maximum number of valleys (namely ten in our code) has been reached; (iii) increasing the number of TAP-like states would lead to too similar valleys, that is some of the overlaps of the TAP-like states would become greater than $`q^{}`$.

At the end of the procedure, the algorithm gives us valleys (families of configurations) and their associated TAP-like states. Not surprisingly, for instances where the clustering leads to several sizeable valleys, the $`P_J(q)`$ has several peaks, and the reverse is also true. Furthermore, we checked the property previously mentioned in section III, namely that configurations in valleys with small mutual overlaps usually differ by sponge-like clusters.

To give the reader a bit of intuition about how the algorithm behaves, let us show what happens for the disorder instance used in figure 1. For that instance, there are two peaks in $`P(q)`$, one near $`q=0.6`$ and the other near $`q=0.8`$. Thus we expect to find two valleys with self-overlap close to $`0.8`$ and whose mutual overlap should be near $`0.6`$. Recall that $`q^{}`$ acts as a resolution; let us then follow the result of the clustering as we go from $`q^{}=0`$ (low resolution) to $`q^{}=1`$ (high resolution).

For small values of $`q^{}`$, the algorithm puts nearly all the configurations into one huge valley. (Recall that we work modulo the $`\pm 1`$ symmetry; the Hamiltonian has no magnetic field and our clustering procedure respects that symmetry by always working modulo this symmetry. In effect, each configuration can be identified with a pair of opposite configurations.) The self-overlap of the corresponding TAP-like state is quite stable, going from $`0.741`$ at $`q^{}=0`$ to $`0.752`$ at $`q^{}=0.61`$. Then at $`q^{}=0.62`$, this big valley breaks up into two sub-valleys. The break-up is accompanied by a discontinuity in the self-overlap: for the largest valley, it jumps to $`0.86`$. Also, the overlap between the two TAP-like states is $`0.61`$. This value is no surprise as the algorithm forbids TAP-like states with overlaps greater than $`q^{}`$; when new valleys appear, one of the inter-valley overlaps should correspond precisely to the cut-off.

Now we increase $`q^{}`$ further. For the range $`0.62q^{}0.72`$, the assignment of the $`2200`$ configurations to the $`2`$ valleys is $`q^{}`$-independent. Beyond $`q^{}=0.72`$, things are qualitatively similar except that a small number of configurations leave these valleys. This can be understood quite simply: the configurations in the valley that are the furthest from the TAP-like center will be pushed out first. Finally, when $`q^{}>0.82`$, the largest valley breaks up and we end up with more small valleys. Clearly, once $`q^{}`$ is close to $`q_{EA}`$, one will not keep a small number of valleys; such values of $`q^{}`$ should not be used. Even in the case of systems without frustration and disorder, this would lead to multiple valleys.

## V Valley-to-valley energy fluctuations

The main goal of this work is to understand how valleys differ; to this end it is convenient, among the $`512`$ disorder instances we have produced, to focus only on those that give rise to more than one valley. For each instance of that type, the $`2200`$ equilibrium configurations cluster into two or more valleys as generated by our algorithm. (For the data shown hereafter, the resolution is $`q^{}=0.5`$.) However, sometimes one of these valleys will be of very small size, containing $`1`$ or $`2`$ configurations. Since we want to investigate typical properties, we have decided to study only valleys of an acceptable size. Quantitatively, we do this by restricting ourselves to the two largest valleys and by demanding that each of them include a substantial fraction of the $`2200`$ spin configurations.

The analysis we present in this section has been performed using those disorder instances where the two largest valleys together contain at least $`80\%`$ of the configurations, and the second largest valley at least $`20\%`$. Out of the $`512`$ disorder instances, these selection criteria left us with $`81`$ instances for $`L=6`$, $`91`$ instances for $`L=8`$, and $`91`$ instances for $`L=12`$. We have repeated the analysis with different thresholds and have found that our conclusions are unaffected by the change. We now begin by discussing the thermodynamic properties of these valleys, while in the sections thereafter we will focus on the corresponding TAP-like states.

The free-energy $`F_r`$ of a valley $`r`$ is defined as

$$e^{F_r/T}={}_{}{}^{}e_{}^{H(C)/T},$$
(4)
where $`^{}`$ indicates that only those configurations $`C`$ belonging to valley $`r`$ are to be included. Note that this sum is the weight of that valley, equal to the valley’s contribution to the partition function $`Z_J`$. Since the Monte Carlo samples configurations according to their Boltzmann factor, it visits the different valleys according to their contribution to $`Z_J`$. Thus, the number of times the Monte Carlo generates a configuration in valley $`r`$ is proportional to $`\mathrm{exp}(F_r/T)`$.

For convenience, we order our valleys according to their size (the largest valley includes $`|V_1|`$ configurations, the next largest $`|V_2|`$, and so on). As a result, the valleys are ordered by increasing free-energies. Defining $`\mathrm{\Delta }F=F_2F_10`$, we can estimate this difference from our data as

$$\mathrm{exp}\left[\frac{\mathrm{\Delta }F}{T}\right]=\frac{|V_2|}{|V_1|}.$$
(5)
In the mean field picture, $`\mathrm{\Delta }F=O(1)`$ with a finite probability at large $`L`$. It is not difficult to see that our selection criteria impose $`\mathrm{\Delta }F=O(1)`$; what is then relevant is the fraction of disorder instances that pass these selection criteria. These fractions are $`15\%`$, $`17\%`$, and $`17\%`$ for $`L=6`$, $`8`$, and $`12`$. Within the error bars, these fractions are stable, in agreement with the mean field picture.

Another prediction of the mean field theory is that $`\mathrm{\Delta }F`$ has an exponential distribution, with a slope becoming small as $`q^{}`$ decreases. In figure 4 we show the numerical values we have found for the probability distribution of $`\mathrm{\Delta }F`$ at $`L=6`$, $`8`$, and $`12`$. The distribution seems quite insensitive to $`L`$ and its shape is rather flat. Again we interpret the data as being compatible with the mean field picture.

Each valley can be thought of as a finite volume equilibrium state; intensive observables are expected to have the same values across different valleys, but extensive quantities can very well fluctuate significantly. Let us thus consider the (extensive) internal energy of valley $`r`$: $`E_r`$ is defined as the average of the energies of the configurations belonging to that valley. We have measured the statistical properties of $`\mathrm{\Delta }EE_2E_1`$. First, we find that the mean of this random variable is positive; that can be understood by saying that since by construction $`\mathrm{\Delta }F>0`$, all other things being equal, one expects $`\mathrm{\Delta }E>0`$.

Second, the spread of the distribution of $`\mathrm{\Delta }E`$ grows with $`L`$; this is very visible in figure 5 and is to be contrasted with the behavior of $`\mathrm{\Delta }F`$. To make this point more quantitative, we have calculated the variances of $`\mathrm{\Delta }E`$; these are given in table II, confirming that the distribution of $`\mathrm{\Delta }E`$ broadens as $`L`$ increases. Third, we find that the linear correlation coefficient of $`\mathrm{\Delta }E`$ and $`\mathrm{\Delta }F`$ is small. This coefficient is defined as

$$C(\mathrm{\Delta }E,\mathrm{\Delta }F)=\frac{\overline{(\mathrm{\Delta }E\overline{\mathrm{\Delta }E})(\mathrm{\Delta }F\overline{\mathrm{\Delta }F})}}{\sigma (\mathrm{\Delta }E)\sigma (\mathrm{\Delta }F)},$$
(6)
where $`\sigma `$ is the standard deviation of the variable of interest and an over-bar stands for the disorder average. The values of $`C`$ are also given in the table. It would be nice to test the possibility that the variance of $`\mathrm{\Delta }E`$ grows as a power of $`L`$, but our range of lattice sizes is too small to seriously test this hypothesis. Nevertheless, the picture we extract from our analysis is that these “large” fluctuations in energy must be compensated by similar fluctuations in entropy, otherwise one could not have $`\mathrm{\Delta }F=O(1)`$. Most probably, these fluctuations play a crucial role in all frustrated disordered systems and they deserve further investigation.

## VI Inter-valley differences: spin and link overlaps

We focus here on the relation between spin and link overlaps: are these two fluctuating quantities correlated? A motivation for this question comes from mean field theory where $`q_l`$ becomes a deterministic function of $`q`$ in the thermodynamic limit. To investigate this possibility, we look at the variance of $`q_l`$ at fixed $`q`$. There are several approaches: the first considers all fluctuations; the second focuses on instance to instance fluctuations only; the third is based on the fluctuations given by the TAP-like states of our valleys.

Consider the first approach. We take equilibrium configurations at random, and determine (for $`q`$ in a small range) the distribution of $`q_l`$ for $`L=6`$, $`8`$, and $`12`$. This distribution is averaged over the disorder, and we consider the resulting variance. We find that the variance decreases as one increases $`L`$; for example, in the bin $`0q0.1`$, the standard deviation is $`0.043`$ at $`L=6`$, $`0.033`$ at $`L=8`$, and $`0.028`$ at $`L=12`$. (The other bins give very similar results.) This decrease is unambiguous, and is highly suggestive of an asymptotic decrease to $`0`$ as $`L\mathrm{}`$. Previous work  on ground states lead to this same conclusion.

Fluctuations in $`q_l`$ come from thermal noise within a given instance and from instance to instance fluctuations. In a naïve way, assume that the variances add linearly, $`\sigma ^2(q_l)\sigma _T^2(q_l)+\sigma _J^2(q_l)`$ where the first one is the intra-instance (purely thermal) variance and the second one is the inter-instance variance. The thermal fluctuations are due to thermally excited droplets inside the valleys; these fluctuations are expected to average out in any scenario, so that $`\sigma _T^2(q_l)0`$ at large $`L`$. The measurements of the previous paragraph show that $`\sigma ^2(q_l)`$ decreases as $`L`$ grows, but this can happen even if $`\sigma _J^2(q_l)`$ saturates at a positive value. It is thus appropriate to determine the $`L`$ dependence of $`\sigma _J^2(q_l)`$ by itself. This leads us to our second method whereby we find the mean of $`q_l`$ in a given disorder instance (and for $`q`$ in a given bin) and then we consider its instance to instance fluctuations. Naturally, the variances now are smaller than in the previous measurements. Is there a sign of saturation? Let us again give our results for $`0q0.1`$: at $`L=6`$, $`\sigma _J(q_l)=0.034`$, at $`L=8`$, $`\sigma _J(q_l)=0.028`$, and at $`L=12`$, $`\sigma _J(q_l)=0.022`$. We have the same decreasing trend as before, strengthening the claim that $`q`$ and $`q_l`$ are related deterministically at large $`L`$.

In the second method just described, we removed the thermal noise in $`q_l`$; can one also remove the thermal noise in $`q`$? Clearly this can be meaningful only if we restrict the averaging so that one stays in given valleys. A simple approach is then to take the mean of $`q`$ for all pairs of configurations belonging to two different but fixed valleys. It is not difficult to see that the corresponding mean overlap is equal to the overlap of the TAP-like states of the two valleys. Let these states be $`^{(1)}`$ and $`^{(2)}`$, having magnetizations $`\{m_i^{(1)}\}`$ and $`\{m_i^{(2)}\}`$; the valley-to-valley spin overlaps average to

$$q(^{(1)},^{(2)})\frac{1}{N}\underset{i}{}m_i^{(1)}m_i^{(2)}$$
(7)
exactly. One is then tempted to extend the analysis by considering the link overlap of these states

$$q_l(^{(1)},^{(2)})\frac{_{ij}m_i^{(1)}m_i^{(2)}m_j^{(1)}m_j^{(2)}}{_{ij}1}$$
(8)
and seeing whether these two overlaps are related deterministically at large $`L`$. We believe that such a test would be appropriate for systems with one-step RSB. But here we have continuous RSB so our TAP-like states are associated with valleys that can have sub-valleys; then, even if there is a deterministic relation between $`q`$ and $`q_l`$ at the level of individual configurations, that will no longer be the case when considering $`q(^{(1)},^{(2)})`$ and $`q_l(^{(1)},^{(2)})`$. Indeed, the averaging over the sub-valleys introduces intrinsic fluctuations that will not go to $`0`$ as $`L\mathrm{}`$. So instead, let us look at the problem differently. To remove the thermal fluctuations, we should remove the droplet excitations. Suppose we think of a valley as a “reference configuration” dressed by any number of droplet excitations. Granted, this picture is simple-minded, but it provides a useful framework. When computing the TAP-like state of that valley, we start with the reference configuration and let the thermal fluctuations decrease its magnetizations from $`\pm 1`$. This suggests that if we bring each $`m_i`$ back to its original value, we will reconstruct this reference configuration which has no thermal noise in it. Our procedure for approximating this reference configuration is simply to project each $`m_i`$ of the TAP-like state back to $`\pm 1`$ according to its sign; we call the resulting configuration the “projected TAP-like state”. It plays the role of a reference configuration, a kind of “heart” of the valley.

In figure 6 we show the collection of spin and link overlaps among these projected TAP-like states in our data sample. The data are scattered, but there is also a clear trend. In figure 7 we show the averaged data, restricted to $`L=12`$, with binning in $`q`$. Superposed in this figure is the best quadratic fit to the data (a quadratic dependence is what one has in the mean field theory). The fit is very satisfactory.

Two questions of interest are:

1. What is the limiting shape of this curve at large $`L`$?
2. In fig. 6, does the scatter go to $`0`$ when $`L`$ grows?

For the limiting curve, we expect it to be a monotonic function, $`q_l`$ growing with $`q`$. For conciseness, let us then just give the values for $`q_l`$ at $`q=0`$. To estimate this mean, we consider all events in the window $`0q0.1`$, leading to $`q_l`$ equal to $`0.428\pm 0.007`$ at $`L=6`$, $`0.483\pm 0.005`$ at $`L=8`$, and $`0.555\pm 0.004`$ at $`L=12`$. The trend is towards increasing values with $`L`$, just as was found for low-energy excitations above the ground state  and from Monte Carlo simulations . Note also that the numbers themselves are very close to those obtained by Houdayer et al. .

To address the second question, we have measured the variance of $`q_l`$ when $`q`$ is restricted to a given range. As above, consider the window $`0q0.1`$; in this range, the standard deviation of $`q_l`$ is $`0.045`$ at $`L=6`$, $`0.037`$ at $`L=8`$, and $`0.031`$ at $`L=12`$. These variances decrease with $`L`$. Again, our range in $`L`$ is too small for us to meaningfully fit these values to a constant plus power law in $`1/L`$. The main point is that the trend is consistent with what was found previously at zero and finite temperature, and gives one confidence that in the large $`L`$ limit $`q`$ and $`q_l`$ are related by a deterministic law.

## VII Inter-valley differences: position space picture

In our working conditions, i.e., at quite low $`T`$, we find that many of the magnetizations $`m_i`$ defining the TAP-like states are close to $`\pm 1`$. This suggests that the projection to $`\pm 1`$ is a reasonable method to reach a valley’s heart. Furthermore, with this projection, we can compare with previous results at $`T=0`$. At stake are the spatial properties of the differences of projected TAP-like states.

To motivate our measurements, let us first “see” what these differences actually look like. In figure 8 we show an example of a difference (called $`G`$) at $`L=12`$ for an instance that passed the selection tests described previously.

The cluster is very space filling; to make this clear, we recommend looking at the display on a computer screen while it is being drawn (this gives far more information than the final printed picture where most of the cluster is hidden). Furthermore, the cluster is topologically non-trivial: one can guess at the numerous handles in its inside, and not surprisingly in our analysis it turns out to be sponge-like (it winds around the lattice in the three directions). One of the central questions for these types of objects concerns the scale beyond which (if any) they become homogeneous. This is a difficult question so it is wisest to focus on the cases for which $`q`$ is small: if any kind of homogeneity can be detected, it should be testable there.

As a probe of heterogeneities, we follow Houdayer et al.  and ask whether it is possible to push a needle all the way through the lattice without hitting the boundary of $`G`$ (the difference of the two projected TAP-like states). This has a direct and intuitive interpretation: let each site of $`G`$ be associated with a unit cube that is opaque, as displayed in the figure. We send light onto one face of the cube and ask whether $`G`$ stops all the light or whether some light is able to pass through the cube without being absorbed. If $`G`$ is homogeneous on all scales larger than the lattice spacing, the probability of absorbing all the light goes to $`1`$ as $`L\mathrm{}`$. Because of the up-down symmetry, we ask in fact whether both $`G`$ and its complement absorb all the light. We find, as in the $`T=0`$ landscape studies , that there is almost always a straight line parallel to a given axis that stays entirely in $`G`$ or entirely in its complement. The system is thus not opaque. We can extend this kind of study by generalizing the needle to a pencil or “tube”: now we try to draw not a single line, but $`4`$ parallel lines whose cross-section is a plaquette of our lattice (the basic square formed by $`4`$ bonds connecting nearest neighbors). The fatter the tube, the less likely it will be possible to slide it through without hitting the surface of $`G`$. To compare with the $`T=0`$ study, we have used the same cross-sections as there: tube “2” has the plaquette cross-section with $`4`$ sites; tube “3” has a diamond cross-section with $`5`$ sites, i.e., a central site and its $`4`$ nearest neighbors; finally tube “4” has for its cross section two adjacent plaquettes with a total of $`6`$ sites (the reader can find pictures of these cross-sections in Houdayer et al. .)

Table III gives the fraction of events for which there is at least one tube along the $`x`$ axis which “spans” the whole cube (i.e., is of length $`L`$) without touching the surface of $`G`$ (the tube sites are thus entirely in $`G`$ or in its complement). The data are for those disorder instances with at least two valleys and with the extra constraint that $`0q0.1`$ where $`q`$ is the spin overlap of the projected TAP-like states.

These data are a bit too noisy to extract a reliable trend. At best, we can say that the data are compatible with increasing values of the spanning probabilities as $`L`$ grows; that is the behavior arising in the energy landscape study. This would mean that the difference of two projected TAP-like states is geometrically similar to the difference of two low-lying system-size excited states. These findings should be considered together with the evidence that the clusters are space-filling (see ): a $`2D`$ plane will intersect the clusters with very high probability. Thus on scales comparable to $`L`$, the clusters are space filling, in agreement with the Parisi mean field theory.

## VIII Discussion and conclusions

We believe that the analysis presented in this note sheds significant light on the issue of the behavior of finite dimensional, realistic spin glasses. On the one hand we stress the physical relevance of the picture of space filling sponges; on the other hand our numerical results from low $`T`$ Monte Carlo simulations give evidence for RSB with a mean field like scenario. Perhaps our most important point is that the results of low $`T`$ numerical simulations give the same indications as the ground state computations. This fact attenuates the long-standing doubt: “Are spin glass numerical simulations representative of low $`T`$ physics, or do they just contain artifacts due to the $`T_c`$ critical point?” Our answer is then that indeed recent numerical simulations do give information about the low $`T`$ physics. Similarly, ground state computations then characterize the low $`T`$ physics behavior of these systems.

A second important issue we address concerns entropic effects. We find that the mean field like features of the system arise with entropy and energy fluctuations; the mutual cancellation of these fluctuations allows the near degeneracy in the free-energy of different states. To see these fluctuations, it is necessary to introduce valleys that are in a sense the building blocks of what can be called finite volume equilibrium states.

A third point worth mentioning is that definitions such as sponge-like or droplet-like differences are very useful; not only do they seem more and more clear-cut as $`L`$ grows, but they also allow for a finer analysis. An illustration of this potential was given via a modified spin overlap (restricted to sponges) that appears to be more robust than the usual one. For the future, we expect these definitions to be extremely useful when trying to understand putative deviations of realistic spin glasses from the standard mean field picture.

There is a lot more that can be done; many issues should be understood in greater detail and checked on larger lattices (our largest size is $`L=12`$). Let us just give a bird’s-eye view of what we feel are some of the main questions. (1): All of our work was for a given temperature; the effect of raising and lowering the temperature should be investigated. (2): The study of the structure of valleys has only begun. Our focus here was on large and well separated valleys; clearly a complete picture requires one to consider the more general case. (3): We have analyzed the EA model with binary ($`\pm `$) couplings. For the lattice sizes we used, this discreteness showed up on the mean energy differences of valleys (see figure 5). Does this effect go away for large sizes and can one confirm numerically the expectation that the Gaussian couplings lead to similar results at large $`L`$?

Acknowledgments

We thank A. Billoire, J.-P. Bouchaud, E. Domany, J. Houdayer, F. Krzakala, M. Mézard and G. Parisi for very stimulating discussions. The LPTMS is an Unité de Recherche de l’Université Paris XI associée au CNRS. One of us (EM) warmly acknowledges the hospitality of the LPTMS during which part of this work was done. The numerical simulations were run on Kalix2, a Linux cluster located at the Physics Department of Cagliari University, and funded by MURST COFIN 1998 (Italy). FZ is supported by an EEC Marie Curie fellowship (contract HPMFCT-2000-00553).

Appendix

We want to characterize valleys or families of spin configurations by their local magnetizations, “à la TAP”, in analogy with what occurs in mean field models. However, from a thermodynamical point of view, it is not easy to define TAP states in a finite volume. Because of this, our definition of valleys is obtained by clustering according to the spin overlaps. One can consider the formal approach of clustering over all of configuration space where each configuration has a given probability, or one can use a more “hands-on” description, where one clusters equilibrium configurations (as obtained from our simulation). Since it is simpler to explain our clustering in the latter framework, we limit ourselves to that case here.

To cluster our configurations, we will use the mutual overlaps and a cut-off parameter $`q^{}`$. To begin, we give two definitions. A valley is just a set of configurations; these will be selected by the clustering algorithm, but what matters is that for each valley $`k`$ one has a collection of configurations $`𝒞_n^k`$ ($`n=1,\mathrm{}|V_k|`$). Second, the TAP-like state associated with valley $`k`$ is the list of magnetizations $`m_i^{(k)}[1,1]`$, one for each site $`i`$ of the lattice. By definition, we take $`m_i^{(k)}`$ to be the mean of $`S_i`$ when considering all the configurations belonging to valley $`k`$.

The logic of our approach can be summarized as follows: (i) a valley defines its TAP-like state; (ii) the TAP-like states tell us whether or not a new configuration should be assigned to one of the $`k`$ valleys. A sensible assignment rule is to impose that the overlap of this new configuration with the corresponding TAP-like state should be sufficiently large. In our algorithm, we first compute the spin overlaps between the new configuration and the $`k`$ TAP-like states; let $`q^{(1)}`$, $`q^{(2)}`$, … $`q^{(k)}`$ be these overlaps. The new configuration is then assigned to valley $`j`$ if $`q^{(j)}`$ is the largest of these overlaps and if $`q^{(j)}q^{}`$.

With these definitions, we can now explain our clustering procedure which is based on defining valleys and their TAP-like states self-consistently. Suppose we have built $`k`$ valleys so far; let their TAP-like states be $`^{(1)}`$, $`^{(2)}`$, …, $`^{(k)}`$. Suppose we want to increase $`k`$ by $`1`$. First, we initialize $`^{(k+1)}`$ by setting its $`m_i^{(k+1)}`$ to be equal to spin values $`S_i`$ of one of the configurations that is not assigned to any valley. Then we go through the rest of the un-assigned configurations, adding them successively to the pool defining the ($`k+1`$)th valley if their overlap $`q^{(k+1)}`$ with $`^{(k+1)}`$ is greater than $`q^{}`$ (the overlap is defined in the usual way, $`q=_1^NS_im_i/N`$). Every time a configuration is added to the pool, we update $`^{(k+1)}`$ so that it always gives the average magnetization at each site for the configurations in valley $`k+1`$.

At the end of this pass, we have $`^{(1)}`$, $`^{(2)}`$, …, $`^{(k+1)}`$. But because $`^{(k+1)}`$ changes as one adds configurations to the pool, some of such configurations may see their overlap $`q^{(k+1)}`$ go below $`q^{}`$. Furthermore, some of the configurations assigned to valleys $`j`$ ($`jk`$) may end up being closer to $`^{(k+1)}`$ than to $`^{(j)}`$. It is thus necessary to rebuild all the valleys, so we first de-assign all the configurations and then we perform a new pass where a configuration is assigned to valley $`r`$ if both $`q^{(r)}>q^{}`$ and $`q^{(r)}>q^{(s)}`$ for all $`sr`$. After three iterations of this build/un-build process, our algorithm takes these to be the valleys and computes the corresponding ($`k+1`$) TAP-like states. These additional iterations improve the self-consistency and reduce the valleys’ dependencies on the ordering of the configurations. This procedure is very close in spirit to that given by Iba and Hukushima ; they also have TAP-like states and valleys that are defined self-consistently by an iterative algorithm.

The last issue pertinent to our method concerns the total number $`M`$ of valleys to be generated. In our implementation, we continue creating more valleys until either of the three following criteria are met: (1) all configurations are assigned to valleys; (2) one has reached the maximum number of valleys (set to $`10`$ in our code); (3) adding one more valley and its TAP-like state would lead to ($`M+1`$) TAP-like states whose cross-overlaps would not all be less than $`q^{}`$. When that happens, we consider that a valley has been subdivided too far and we do not accept this ($`M+1`$)th TAP-like state.