format: text

text:
# The dynamics of quantum phases in a spinor condensate

## I Introduction

Since the observation of Bose-Einstein condensation in trapped dilute alkali atoms bec ; mit ; rice , the coherence properties of the condensate has become the focus of many theoretical studies lewenstein ; wright ; imamoglu ; villain97 ; law ; villain99 . As was initially pointed out in Refs. lewenstein ; wright , the number fluctuations associated with the assumption of a coherent state macroscopic condensate wavefunction leads to the “diffusion” (or spreading) of the initial phase. For a scalar, one component atomic condensate, this diffusion, the attempt to restore the $`U(1)`$ symmetry of the underline interacting system, can be formulated in terms of the dynamics of a zero mode, or the Goldstone mode lewenstein . More physically meaningful discussions in terms of the relative phases of two condensates were studied soon afterwards imamoglu ; villain97 ; law ; villain99 . Starting with the remarkable direct observation of the first order coherence in an interference experiment mitint , direct relationship between the number fluctuation and quantum phase dynamics was first observed with a condensate in a periodic potential kasevich , and more recently, in the remarkable Mott insulating state hansch ; zoller .

The emergence of spinor condensates kurn ; mike  (of atoms with hyperfine quantum number $`F=1`$) has created new opportunities to understand quantum coherence and the associated phase dynamics within a three component condensate burnett ; ashab . In an earlier paper yil , we have presented important physical insight into quantum phase diffusions of a spinor condensate. We have investigated the coupled zero mode dynamics as a consequence of the usual mean field treatment that calls for the breaking of two continuous symmetries: a U(1) gauge transformation $`e^{i\theta }`$ and SO(3) spin rotations $`𝒰(\alpha ,\beta ,\tau )=e^{iF_z\alpha }e^{iF_y\beta }e^{iF_z\tau }`$ (in the absence of an external magnetic field) ho .

This paper provides a detailed investigation of the quantum phase dynamics for a spinor-1 condensate. As before, we will focus on the case of ferromagnetic interactions, for which the validity of the single spatial mode approximation proves to be a convenient starting point for the emergence of transparent physical illustrations of the zero mode dynamics su . The fluctuations of atom numbers in each of the spin states will be shown to be connected directly, in this case, to the fluctuations of the pointing direction of the macroscopic spin (or magnetic dipole moment) of the condensate. This paper is organized as follows. In section II, after introducing the Hamiltonian for a spinor-1 condensate, we explicitly work out the Heisenberg equations for the atomic field operators. Following the standard Bogoliubov approach bog , we obtain the coupled Gross-Pitaevskii equations for the mean field of the condensate and the coupled Bogoliubov-de Gennes equations for the quantum fluctuations. We then generalize the Hermitian operators for the condensate number and phase fluctuations (of the Goldstone mode) as introduced earlier lewenstein  to each of the spin component, and derive their dynamic equations under the rotating wave function approximation. In section III, we perform a detailed study of the stationary phase diffusion of a spinor-1 condensate within the single mode approximation valid exactly for ferromagnetic interactions. Section IV addresses the influence of an external magnetic field, and in section V we discuss the same phase dynamics for a condensate with anti-ferromagnetic interactions. We conclude in section VI.

## II Formulation

We consider a system of $`N`$ spin-1 bosonic atoms interacting via only s-wave scattering ho ; ohmi ; goldstein ; pu . In the second-quantized form, the Hamiltonian of our system becomes ho ; ohmi

$`H`$ $`=`$ $`{\displaystyle \underset{i}{}}{\displaystyle 𝑑\stackrel{}{r}\psi _i^{}(\stackrel{}{r})\left[\frac{\mathrm{}^2^2}{2M}+V_{\mathrm{ext}}(\stackrel{}{r})\right]\psi _i(\stackrel{}{r})}`$ (1)
$`+{\displaystyle \frac{c_0}{2}}{\displaystyle \underset{i,j}{}}{\displaystyle 𝑑\stackrel{}{r}\psi _i^{}(\stackrel{}{r})\psi _j^{}(\stackrel{}{r})\psi _i(\stackrel{}{r})\psi _j(\stackrel{}{r})}`$
$`+{\displaystyle \frac{c_2}{2}}{\displaystyle \underset{i,j,k,l}{}}{\displaystyle 𝑑\stackrel{}{r}\psi _i^{}(\stackrel{}{r})\psi _j^{}(\stackrel{}{r})\stackrel{}{𝐅}_{ik}\stackrel{}{𝐅}_{jl}\psi _l(\stackrel{}{r})\psi _k(\stackrel{}{r})},`$
where $`\psi _j(\stackrel{}{r})`$ ($`j=+,0,`$) denotes the annihilation operator for the $`j`$-th component of a spinor-1 field. The external trapping potential $`V_{\mathrm{ext}}(\stackrel{}{r})`$ is assumed spin-independent as in a far off-resonant optical dipole force trap (FORT) so that atomic spin degrees of freedom are completely accessible. The pair interaction coefficients are $`c_0=4\pi \mathrm{}^2(a_0+2a_2)/3M`$ and $`c_2=4\pi \mathrm{}^2(a_2a_0)/3M`$, with $`a_0`$ ($`a_2`$) the s-wave scattering length for two spin-1 atoms in the combined symmetric channel of total spin $`0`$ ($`2`$), and $`M`$ is the mass of atom. $`\stackrel{}{𝐅}`$ is spin 1 matrix representation. The interacting terms can be regrouped to show that they include respectively self-scattering, cross-scattering, and the spin-relaxation pu .

From the Hamiltonian (1), one can easily derive the Heisenberg equations for the field operators

$`i\mathrm{}\dot{\psi }_+`$ $`=`$ $`\left[{\displaystyle \frac{\mathrm{}^2^2}{2M}}+V_{\mathrm{ext}}(\stackrel{}{r})\right]\psi _++c_2\psi _{}^{}\psi _0^2`$
$`+c_0(\psi _+^{}\psi _++\psi _0^{}\psi _0+\psi _{}^{}\psi _{})\psi _+`$
$`+c_2(\psi _+^{}\psi _++\psi _0^{}\psi _0\psi _{}^{}\psi _{})\psi _+,`$
$`i\mathrm{}\dot{\psi }_0`$ $`=`$ $`\left[{\displaystyle \frac{\mathrm{}^2^2}{2M}}+V_{\mathrm{ext}}(\stackrel{}{r})\right]\psi _0+2c_2\psi _0^{}\psi _+\psi _{}`$
$`+c_0(\psi _+^{}\psi _++\psi _0^{}\psi _0+\psi _{}^{}\psi _{})\psi _0`$
$`+c_2(\psi _+^{}\psi _++\psi _{}^{}\psi _{})\psi _0,`$
$`i\mathrm{}\dot{\psi }_{}`$ $`=`$ $`\left[{\displaystyle \frac{\mathrm{}^2^2}{2M}}+V_{\mathrm{ext}}(\stackrel{}{r})\right]\psi _{}+c_2\psi _+^{}\psi _0^2`$ (2)
$`+c_0(\psi _+^{}\psi _++\psi _0^{}\psi _0+\psi _{}^{}\psi _{})\psi _{}`$
$`+c_2(\psi _{}^{}\psi _{}+\psi _0^{}\psi _0\psi _+^{}\psi _+)\psi _{}.`$
Following the Hartree mean-field theory we assume there are three ‘large’ condensate components $`\varphi _j`$ around which we can study the small quantum fluctuations (off-condensate excitations) via the Bogoluibov approximation

$`\psi _j(\stackrel{}{r})=\sqrt{N}\varphi _j(\stackrel{}{r})+\delta \psi _j(\stackrel{}{r}),`$ (3)
with $`N`$ the total number of atoms. We note that such an approach can become questionable in a multi-component system where certain component becomes negligibly small and its quantum nature becomes important. For ferromagnetic interactions, however, our previous studies su  show that the ground state wave functions can always be expressed as $`\varphi _j\varphi `$ and hence the Bogoluibov prescription can always be applied. By substituting Eq. (3) into (2), we find in the zeroth order, the coupled Gross-Pitaevskii equations (GPEs),

$`i\mathrm{}\dot{\varphi }_+`$ $`=`$ $`\left[+c_0\rho +c_2(\rho _++\rho _0\rho _{})\right]\varphi _++c_2N\varphi _{}^{}\varphi _0^2,`$
$`i\mathrm{}\dot{\varphi }_0`$ $`=`$ $`\left[+c_0\rho +c_2(\rho _++\rho _{})\right]\varphi _0+2c_2N\varphi _0^{}\varphi _+\varphi _{},`$
$`i\mathrm{}\dot{\varphi }_{}`$ $`=`$ $`\left[+c_0\rho +c_2(\rho _{}+\rho _0\rho _+)\right]\varphi _{}+c_2N\varphi _+^{}\varphi _0^2,`$
with $`=\mathrm{}^2^2/2M+V_{\mathrm{ext}}(\stackrel{}{r})`$, the $`j`$-th component condensate density $`\rho _j(\stackrel{}{r})=N|\varphi _j(\stackrel{}{r})|^2`$, and the total condensate density $`\rho (\stackrel{}{r})=\rho _+(\stackrel{}{r})+\rho _0(\stackrel{}{r})+\rho _{}(\stackrel{}{r})`$. $`𝑑\stackrel{}{r}\rho _j(\stackrel{}{r})=N_j`$ is the number of condensed atoms in the $`j`$-th component. In the first order, the quantum fluctuations obey the following set of Bogoluibov-de Gennes equations (BdGEs),

$`i\mathrm{}{\displaystyle \frac{\delta \psi _+}{t}}`$ $`=`$ $`\left[+c_0(\rho +\rho _+)+c_2(2\rho _++\rho _0\rho _{})\right]\delta \psi _+`$
$`+N(c_+\varphi _0^{}\varphi _++2c_2\varphi _{}^{}\varphi _0)\delta \psi _0`$
$`+c_{}N\varphi _{}^{}\varphi _+\delta \psi _{}+c_+N\varphi _+^2\delta \psi _+^{}`$
$`+c_+N\varphi _+\varphi _0\delta \psi _0^{}+N(c_{}\varphi _+\varphi _{}+c_2\varphi _0^2)\delta \psi _{}^{},`$
$`i\mathrm{}{\displaystyle \frac{\delta \psi _0}{t}}`$ $`=`$ $`\left[+c_0(\rho +\rho _0)+c_2(\rho _++\rho _{})\right]\delta \psi _0`$
$`+N(c_+\varphi _+^{}\varphi _0+2c_2\varphi _0^{}\varphi _{})\delta \psi _+`$
$`+N(c_+\varphi _{}^{}\varphi _0+2c_2\varphi _0^{}\varphi _+)\delta \psi _{}`$
$`+c_+N\varphi _+\varphi _0\delta \psi _+^{}+c_+N\varphi _0\varphi _{}\delta \psi _{}^{}`$
$`+N(c_0\varphi _0^2+2c_2\varphi _+\varphi _{})\delta \psi _0^{},`$
$`i\mathrm{}{\displaystyle \frac{\delta \psi _{}}{t}}`$ $`=`$ $`\left[+c_0(\rho +\rho _{})+c_2(2\rho _{}+\rho _0\rho _+)\right]\delta \psi _{}`$
$`+N(c_+\varphi _0^{}\varphi _{}+2c_2\varphi _+^{}\varphi _0)\delta \psi _0`$
$`+c_{}N\varphi _+^{}\varphi _{}\delta \psi _++c_+N\varphi _{}^2\delta \psi _{}^{}`$
$`+c_+N\varphi _{}\varphi _0\delta \psi _0^{}+N(c_{}\varphi _+\varphi _{}+c_2\varphi _0^2)\delta \psi _+^{},`$
where we have defined $`c_\pm c_0\pm c_2`$.

To study the quantum phase dynamics, we introduce the number fluctuation operators lewenstein

$`P_j={\displaystyle 𝑑\stackrel{}{r}\left[\varphi _j^{}(\stackrel{}{r})\delta \psi _j(\stackrel{}{r})+\varphi _j(\stackrel{}{r})\delta \psi _j^{}(\stackrel{}{r})\right]}.`$ (6)
Using the GPEs (LABEL:gpe3) and the BdGEs (LABEL:fluc) above, we find the evolution equations for these operators,

$`i\mathrm{}{\displaystyle \frac{P_+}{t}}`$ $`=`$ $`{\displaystyle \frac{i\mathrm{}}{2}}{\displaystyle \frac{P_0}{t}}=i\mathrm{}{\displaystyle \frac{P_{}}{t}}`$ (7)
$`=`$ $`c_2N{\displaystyle }d\stackrel{}{r}(\varphi _{}^{}\varphi _0^2\delta \psi _+^{}+\varphi _+^{}\varphi _0^2\delta \psi _{}^{}`$
$`+2\varphi _+^{}\varphi _{}^{}\varphi _0\delta \psi _0h.c.).`$
Thus both $`P_{\mathrm{tot}}=P_++P_0+P_{}`$ and $`P_+P_{}`$ are constants of motion, which are in fact obvious since the Hamiltonian (1) conserves the total number of atoms and the total magnetization.

We define the conjugate phase operator $`Q_{\mathrm{tot}}=Q_++Q_0+Q_{}`$ according to

$`Q_j=i\mathrm{}{\displaystyle 𝑑\stackrel{}{r}\left[\theta _j^{}(\stackrel{}{r})\delta \psi _j(\stackrel{}{r})\theta _j(\stackrel{}{r})\delta \psi _j^{}(\stackrel{}{r})\right]}.`$ (8)
The canonical quantization condition $`[Q_{\mathrm{tot}},P_{\mathrm{tot}}]=i\mathrm{}`$ is satisfied if the constraint $`J_++J_0+J_{}1`$ is enforced with $`J_j𝑑\stackrel{}{r}\left[\theta _j^{}(\stackrel{}{r})\varphi _j(\stackrel{}{r})+\theta _j(\stackrel{}{r})\varphi _j^{}(\stackrel{}{r})\right]`$. From the defining equation for $`Q_{\mathrm{tot}}`$,

$`{\displaystyle \frac{dQ_{\mathrm{tot}}}{dt}}=N\stackrel{~}{u}P_{\mathrm{tot}},`$ (9)
the Goldstone mode inertial parameter $`\stackrel{~}{u}`$ can be determined. We obtain the dynamic equations for the phase functions as,

$`i\mathrm{}\dot{\theta }_+`$ $`=`$ $`\left[+c_0(\rho +\rho _+)+c_2(2\rho _++\rho _0\rho _{})\right]\theta _+`$
$`+N(c_+\varphi _0^{}\varphi _++2c_2\varphi _{}^{}\varphi _0)\theta _0+c_{}N\varphi _{}^{}\varphi _+\theta _{}`$
$`+c_+N\varphi _+^2\theta _+^{}+c_+N\varphi _+\varphi _0\theta _0^{}`$
$`+N(c_{}\varphi _+\varphi _{}+c_2\varphi _0^2)\theta _{}^{}N\stackrel{~}{u}\varphi _+,`$
$`i\mathrm{}\dot{\theta }_0`$ $`=`$ $`\left[+c_0(\rho +\rho _0)+c_2(\rho _++\rho _{})\right]\theta _0`$
$`+N(c_+\varphi _+^{}\varphi _0+2c_2\varphi _0^{}\varphi _{})\theta _+`$
$`+N(c_+\varphi _{}^{}\varphi _0+2c_2\varphi _0^{}\varphi _+)\theta _{}`$
$`+c_+N\varphi _+\varphi _0\theta _+^{}+c_+N\varphi _0\varphi _{}\theta _{}^{}`$
$`+N(c_0\varphi _0^2+2c_2\varphi _+\varphi _{})\theta _0^{}N\stackrel{~}{u}\varphi _0,`$
$`i\mathrm{}\dot{\theta }_{}`$ $`=`$ $`\left[+c_0(\rho +\rho _{})+c_2(2\rho _{}+\rho _0\rho _+)\right]\theta _{}`$ (10)
$`+N(c_+\varphi _0^{}\varphi _{}+2c_2\varphi _+^{}\varphi _0)\theta _0+c_{}N\varphi _+^{}\varphi _{}\theta _+`$
$`+c_+N\varphi _{}^2\theta _{}^{}+c_+N\varphi _{}\varphi _0\theta _0^{}`$
$`+N(c_{}\varphi _+\varphi _{}+c_2\varphi _0^2)\theta _+^{}N\stackrel{~}{u}\varphi _{}.`$
Combing Eqs. (LABEL:fluc) and (10), the dynamic equations of $`Q_j`$ can be easily obtained to be

$`\dot{Q}_+`$ $`=`$ $`N\stackrel{~}{u}P_++N{\displaystyle }d\stackrel{}{r}\{[(c_+\varphi _+^{}\varphi _0+2c_2\varphi _0^{}\varphi _{})\theta _0^{}+c_+\varphi _+^{}\varphi _0^{}\theta _0+(c_{}\varphi _+^{}\varphi _{}^{}+c_2\varphi _0^2)\theta _{}+c_{}\varphi _+^{}\varphi _{}\theta _{}^{}]\delta \psi _+`$
$`+(c_+\varphi _0^{}\varphi _+\theta _+^{}+c_+\varphi _+^{}\varphi _0^{}\theta _++2c_2\varphi _{}^{}\varphi _0\theta _+^{})\delta \psi _0+(c_{}\varphi _+^{}\varphi _{}^{}\theta _++c_{}\varphi _{}^{}\varphi _+\theta _+^{}+c_2\varphi _0^2\theta _+)\delta \psi _{}+h.c\},`$
$`\dot{Q}_0`$ $`=`$ $`N\stackrel{~}{u}P_0+N{\displaystyle }d\stackrel{}{r}\{[(c_+\varphi _0^{}\varphi _++2c_2\varphi _{}^{}\varphi _0)\theta _+^{}+c_+\varphi _+^{}\varphi _0^{}\theta _++(c_+\varphi _0^{}\varphi _{}+2c_2\varphi _+^{}\varphi _0)\theta _{}^{}+c_+\varphi _0^{}\varphi _{}^{}\theta _{}]\delta \psi _0`$
$`+(c_+\varphi _+^{}\varphi _0\theta _0^{}+c_+\varphi _+^{}\varphi _0^{}\theta _0+2c_2\varphi _0^{}\varphi _{}\theta _0^{})\delta \psi _++(c_+\varphi _{}^{}\varphi _0\theta _0^{}+c_+\varphi _0^{}\varphi _{}^{}\theta _0+2c_2\varphi _0^{}\varphi _+\theta _0^{})\delta \psi _{}+h.c\},`$
$`\dot{Q}_{}`$ $`=`$ $`N\stackrel{~}{u}P_{}+N{\displaystyle }d\stackrel{}{r}\{[(c_+\varphi _{}^{}\varphi _0+2c_2\varphi _0^{}\varphi _+)\theta _0^{}+c_+\varphi _{}^{}\varphi _0^{}\theta _0+(c_{}\varphi _+^{}\varphi _{}^{}+c_2\varphi _0^2)\theta _++c_{}\varphi _{}^{}\varphi _+\theta _+^{}]\delta \psi _{}`$ (11)
$`+(c_+\varphi _0^{}\varphi _{}\theta _{}^{}+c_+\varphi _{}^{}\varphi _0^{}\theta _{}+2c_2\varphi _+^{}\varphi _0\theta _{}^{})\delta \psi _0+(c_{}\varphi _+^{}\varphi _{}^{}\theta _{}+c_{}\varphi _+^{}\varphi _{}\theta _{}^{}+c_2\varphi _0^2\theta _{})\delta \psi _++h.c\}.`$
We note that $`[Q_j,P_k]=i\mathrm{}\delta _{jk}J_j`$, $`[\delta \psi _j(\stackrel{}{r}),P_k]=\delta _{jk}\varphi _j(\stackrel{}{r})`$, and $`[\delta \psi _j(\stackrel{}{r}),Q_k]=i\mathrm{}\delta _{jk}\theta _j(\stackrel{}{r})`$. This prompts us to make the rotating wave approximation (RWA) villain97 ; villain99  by assuming that $`\delta \psi _j(\stackrel{}{r})=_ja_j(\stackrel{}{r})P_j+_jb_j(\stackrel{}{r})Q_j`$, which leads to $`a_j(\stackrel{}{r})=\theta _j(\stackrel{}{r})/J_j`$ and $`b_j(\stackrel{}{r})=\varphi _j(\stackrel{}{r})/i\mathrm{}J_j`$. Therefore, under the RWA, we obtain $`\delta \psi _j(\stackrel{}{r})=\theta _j(\stackrel{}{r})P_j/J_j+\varphi _j(\stackrel{}{r})Q_j/i\mathrm{}J_j`$. With this result, the complete dynamic equations for the number and phase fluctuation operators become

$`\dot{P}_+`$ $`=`$ $`{\displaystyle \frac{\dot{P}_0}{2}}=\dot{P}_{}`$
$`=`$ $`{\displaystyle \frac{2c_2N}{\mathrm{}}}\left[{\displaystyle \frac{(\varphi _{}\theta _+|\varphi _0^2)^{\prime \prime }}{J_+}}P_+2{\displaystyle \frac{(\varphi _0\theta _0|\varphi _+\varphi _{})^{\prime \prime }}{J_0}}P_0+{\displaystyle \frac{(\varphi _+\theta _{}|\varphi _0^2)^{\prime \prime }}{J_{}}}P_{}\right]`$
$`+{\displaystyle \frac{2c_2N}{\mathrm{}^2}}(\varphi _+\varphi _{}|\varphi _0^2)^{}\left({\displaystyle \frac{Q_+}{J_+}}{\displaystyle \frac{2Q_0}{J_0}}+{\displaystyle \frac{Q_{}}{J_{}}}\right),`$
$`\dot{Q}_+`$ $`=`$ $`N\stackrel{~}{u}P_+{\displaystyle \frac{2N}{J_+}}\left[2c_{}I_++2c_+I_{+0}+2c_2(\varphi _0\theta _0|\varphi _{}\theta _+)^{}+c_2(\varphi _0^2|\theta _{}\theta _+)^{}\right]P_+`$
$`+{\displaystyle \frac{4N}{J_0}}\left[c_+I_{+0}+c_2(\varphi _0\theta _0|\varphi _{}\theta _+)^{}\right]P_0+{\displaystyle \frac{2N}{J_{}}}\left[2c_{}I_++c_2(\varphi _0^2|\theta _+\theta _{})^{}\right]P_{}`$
$`+{\displaystyle \frac{2c_2N}{\mathrm{}J_+}}\left[2(\varphi _+\varphi _{}|\varphi _0\theta _0)^{\prime \prime }+(\varphi _+\theta _{}|\varphi _0^2)^{\prime \prime }\right]Q_++{\displaystyle \frac{4c_2N}{\mathrm{}J_0}}(\varphi _{}\theta _+|\varphi _0^2)^{\prime \prime }Q_0+{\displaystyle \frac{2c_2N}{\mathrm{}J_{}}}(\varphi _0^2|\theta _+\varphi _{})^{\prime \prime }Q_{},`$
$`\dot{Q}_0`$ $`=`$ $`N\stackrel{~}{u}P_0{\displaystyle \frac{4N}{J_0}}\left[c_+I_{+0}+c_+I_0+c_2(\varphi _{}\theta _+|\varphi _0\theta _0)^{}+c_2(\varphi _+\theta _{}|\varphi _0\theta _0)^{}\right]P_0`$
$`+{\displaystyle \frac{4N}{J_+}}\left[c_+I_{+0}+c_2(\varphi _0\theta _0|\varphi _{}\theta _+)^{}\right]P_++{\displaystyle \frac{4N}{J_{}}}\left[c_+I_0+c_2(\varphi _0\theta _0|\varphi _+\theta _{})^{}\right]P_{}`$
$`+{\displaystyle \frac{4c_2N}{\mathrm{}J_+}}(\varphi _0\theta _0|\varphi _+\varphi _{})^{\prime \prime }Q_++{\displaystyle \frac{4c_2N}{\mathrm{}J_0}}[(\varphi _0^2|\varphi _{}\theta _+)^{\prime \prime }+(\varphi _0^2|\varphi _+\theta _{})^{\prime \prime }]Q_0+{\displaystyle \frac{4c_2N}{\mathrm{}J_{}}}(\varphi _0\theta _0|\varphi _+\varphi _{})^{\prime \prime }Q_{},`$
$`\dot{Q}_{}`$ $`=`$ $`N\stackrel{~}{u}P_{}{\displaystyle \frac{2N}{J_{}}}\left[2c_{}I_++2c_+I_0+2c_2(\varphi _0\theta _0|\varphi _+\theta _{})^{}+c_2(\varphi _0^2|\theta _+\theta _{})^{}\right]P_{}`$ (12)
$`+{\displaystyle \frac{4N}{J_0}}\left[c_+I_0+c_2(\varphi _0\theta _0|\varphi _+\theta _{})^{}\right]P_0+{\displaystyle \frac{2N}{J_+}}\left[2c_{}I_++c_2(\varphi _0^2|\theta _{}\theta _+)^{}\right]P_+`$
$`+{\displaystyle \frac{2c_2N}{\mathrm{}J_+}}(\varphi _0^2|\theta _{}\varphi _+)^{\prime \prime }Q_++{\displaystyle \frac{4c_2N}{\mathrm{}J_0}}(\varphi _+\theta _{}|\varphi _0^2)^{\prime \prime }Q_0+{\displaystyle \frac{2c_2N}{\mathrm{}J_{}}}\left[2(\varphi _{}\varphi _+|\varphi _0\theta _0)^{\prime \prime }+(\varphi _{}\theta _+|\varphi _0^2)^{\prime \prime }\right]Q_{},`$
where we have introduced the following shorthand notation; $`(fg|hy)=𝑑\stackrel{}{r}f^{}(\stackrel{}{r})g^{}(\stackrel{}{r})h(\stackrel{}{r})y(\stackrel{}{r})`$, $`(.)^{}=\mathrm{Re}(.)`$, $`(.)^{\prime \prime }=\mathrm{Im}(.)`$, and $`I_{ml}=𝑑\stackrel{}{r}\left(\varphi _m\theta _m^{}\right)^{}\left(\varphi _l\theta _l^{}\right)^{}`$. Equation (12) thus completely determines the zero mode condensate fluctuations. Before finding the corresponding compact zero mode Hamiltonian in quadratic forms of $`P_j`$ and $`Q_j`$, we will first attempt to simplify in the next section using general properties of $`\varphi _j(\stackrel{}{r})`$ and $`\theta _j(\stackrel{}{r})`$ for a ferromagnetic spinor-1 condensate su .

## III Phase diffusions and the zero modes

In general, the ground state wave functions of a spinor condensate can be expressed as

$`\varphi _j(\stackrel{}{r},t)=|\varphi _j(\stackrel{}{r})|e^{i\mu t/\mathrm{}+i\alpha _j},`$ (13)
with a common chemical potential $`\mu `$ for all spin components when the external magnetic field is negligible. Based on the minimization for the total energy Eq. (1) within the mean field theory, it was shown in Refs. isoshima ; robins  that there exists an important relation among $`\alpha _j`$, given by $`\alpha _++\alpha _{}2\alpha _0=0`$. We further proved in Ref. su  that the ground state wave function for each of the spin component takes the same spatial shape for (ferromagnetic interactions), i.e.

$`\varphi _j(\stackrel{}{r})=\sqrt{n_j}\varphi (\stackrel{}{r})e^{i\alpha _j},`$ (14)
with the real-valued mode function $`\varphi (\stackrel{}{r})`$ (normalized to unity) governed by an equivalent scalar condensate GPE

$`\left[+c_+N\varphi ^2(\stackrel{}{r})\right]\varphi (\stackrel{}{r})=\mu \varphi (\stackrel{}{r}),`$ (15)
of a scattering length $`a_2`$ (note $`c_+a_2`$). Atoms thus only collide in the symmetric total spin $`F=2`$ channel in a ferromagnetic state in order to maintain their individual spins parallel. $`n_j=N_j/N`$ is the ratio of the number of atoms in the $`i`$-th spin component to the total number of atoms. For any given magnetization $``$, $`n_j`$ is given explicitly by $`n_\pm =(1\pm m)^2/4`$ and $`n_0=(1m^2)/2`$ with $`m=/N`$ pu ; su .

Similarly, we look for phase functions of the form

$`\theta _j(\stackrel{}{r},t)=|\theta _j(\stackrel{}{r})|e^{i\mu _jt/\mathrm{}+i\alpha _j}.`$ (16)
Substituting Eq. (16) into Eq. (10), we find that $`\theta _j(\stackrel{}{r},t)`$ also evolve in time as $`e^{i\mu t/\mathrm{}}`$. To focus on the steady state quantum fluctuation properties, we may therefore neglect the time-dependent part in $`\varphi _j(\stackrel{}{r},t)`$ and $`\theta _j(\stackrel{}{r},t)`$ from now on. The steady state equations for phase functions now become

$`\left[\mu +c_0(\rho +2\rho _+)+c_2(3\rho _++\rho _0\rho _{})\right]\xi _+`$
$`+2N(c_+\varphi _++c_2\varphi _{})\varphi _0\xi _0+N(2c_{}\varphi _+\varphi _{}+c_2\varphi _0^2)\xi _{}`$
$`=N\stackrel{~}{u}\varphi _+,`$
$`\left[\mu +c_0(\rho +2\rho _0)+c_2(\rho _++2N\varphi _+\varphi _{}+\rho _{})\right]\xi _0`$
$`+2N(c_+\varphi _++c_2\varphi _{})\varphi _0\xi _++2N(c_+\varphi _{}+c_2\varphi _+)\varphi _0\xi _{}`$
$`=N\stackrel{~}{u}\varphi _0,`$
$`\left[\mu +c_0(\rho +2\rho _{})+c_2(3\rho _{}+\rho _0\rho _+)\right]\xi _{}`$
$`+2N(c_+\varphi _{}+c_2\varphi _+)\varphi _0\xi _0+N(2c_{}\varphi _+\varphi _{}+c_2\varphi _0^2)\xi _+`$
$`=N\stackrel{~}{u}\varphi _{},`$ (17)
where $`\xi _j(\stackrel{}{r})=\theta _j(\stackrel{}{r})e^{i\alpha _j}`$ are now real-valued functions. One can easily check that similar to the single mode approximation for $`\varphi _j(\stackrel{}{r})`$, $`\xi _j(\stackrel{}{r})`$ can also be expressed as $`\xi _j(\stackrel{}{r})=n_j^{1/2}\theta (\stackrel{}{r})`$ with $`\theta (\stackrel{}{r})`$ governed by the following equation

$`\left[\mu +3c_+N\varphi ^2(\stackrel{}{r})\right]\theta (\stackrel{}{r})=N\stackrel{~}{u}\varphi (\stackrel{}{r}),`$ (18)
and $`\stackrel{~}{u}`$ is determined through the normalization constraint $`𝑑\stackrel{}{r}\varphi (\stackrel{}{r})\theta (\stackrel{}{r})=1/2`$. Therefore, for ferromagnetic interactions, the phase functions can be generally expressed as

$`\theta _j(\stackrel{}{r})=\sqrt{n_j}\theta (\stackrel{}{r})e^{i\alpha _j}.`$ (19)
We note that Eq. (18) is in fact identical to the equation satisfied by the phase function of a scalar condensate lewenstein , a point easily understood for ferromagnetic interactions when all spins align along the same direction. If we take this direction as the quantization axis, then the spinor condensate is essentially a scalar condensate since its total magnetization $``$ is conserved. All overlap integrals in Eq. (12) are presented in Appendix A. They are all real quantities with the choice of phase parameters.

We now attempt to express the zero mode dynamics in terms of its associated Hamiltonian. Since $`P_j`$ and $`Q_j`$ are not canonically conjugated variables, we define $`p_j^{}P_j`$ and $`q_j^{}Q_j/J_j`$ such that they become canonically conjugates. The zero mode Hamiltonian that governs the dynamics of Eq. (12) is

$`{\displaystyle \frac{H_{\mathrm{zero}}}{N}}`$ $`=`$ $`\left({\displaystyle \frac{\stackrel{~}{u}}{2n_+𝒪_{\varphi \theta }}}2c_+{\displaystyle \frac{1n_+}{n_+}}2c_2{\displaystyle \frac{n_{}}{n_+}}\right)p_+^2𝒪_{\varphi \theta }`$ (20)
$`+\left({\displaystyle \frac{\stackrel{~}{u}}{2n_0𝒪_{\varphi \theta }}}2c_+{\displaystyle \frac{1n_0}{n_0}}2c_2\right)p_0^2𝒪_{\varphi \theta }`$
$`+\left({\displaystyle \frac{\stackrel{~}{u}}{2n_{}𝒪_{\varphi \theta }}}2c_+{\displaystyle \frac{1n_{}}{n_{}}}2c_2{\displaystyle \frac{n_+}{n_{}}}\right)p_{}^2𝒪_{\varphi \theta }`$
$`+4\left[c_++c_2\left({\displaystyle \frac{n_{}}{n_+}}\right)^{1/2}\right]p_+^{}p_0^{}𝒪_{\varphi \theta }`$
$`+4\left[c_++c_2\left({\displaystyle \frac{n_+}{n_{}}}\right)^{1/2}\right]p_{}^{}p_0^{}𝒪_{\varphi \theta }`$
$`+4c_0p_+^{}p_{}^{}𝒪_{\varphi \theta }{\displaystyle \frac{c_2n_0^2𝒪_{\varphi \varphi }}{2\mathrm{}^2}}\left(q_+^{}2q_0^{}+q_{}^{}\right)^2`$
$`=`$ $`𝐩^T𝒜𝐩^{}+𝐪^T𝐪^{},`$
where $`𝐩^T=(p_+^{},p_0^{},p_{}^{})`$, $`𝐪^T=(q_+^{},q_0^{},q_{}^{})`$, $`𝒜`$ and $``$ are two Hermitian matrices whose elements are easily identified from above equation. We also present, in Appendix B, the zero mode Hamiltonian in the form of overlap integrals. To simplify notations, we have also defined $`𝒪_{\varphi \theta }𝑑\stackrel{}{r}\varphi ^2(\stackrel{}{r})\theta ^2(\stackrel{}{r})`$ and $`𝒪_{\varphi \varphi }𝑑\stackrel{}{r}\varphi ^4(\stackrel{}{r})`$. We find that matrix $``$ is positive definite since $`c_2<0`$ for ferromagnetic spinor condensates. To show the positive definiteness of $`𝒜`$, we decompose $`\stackrel{~}{u}`$ as

$`\stackrel{~}{u}=\epsilon +4c_+𝒪_{\varphi \theta },`$ (21)
where

$`\epsilon ={\displaystyle \frac{2}{N}}{\displaystyle 𝑑\stackrel{}{r}\theta (\stackrel{}{r})\left[\mu +c_+N\varphi ^2(\stackrel{}{r})\right]\theta (\stackrel{}{r})}.`$ (22)
We note that $`\epsilon 0`$ for noninteracting atoms \[$`\varphi (\stackrel{}{r})=\theta (\stackrel{}{r})`$\]. However, it’s nonzero in general and we can prove that $`\epsilon >0`$ since $`\theta (\stackrel{}{r})`$ is the solution of Eq. (18). By using Eq. (21), we then obtain

$`𝐩^T𝒜𝐩^{}`$ $`=`$ $`{\displaystyle \frac{\epsilon }{2}}\left({\displaystyle \frac{p_+^2}{n_+}}+{\displaystyle \frac{p_0^2}{n_0}}+{\displaystyle \frac{p_{}^2}{n_{}}}\right)+2𝒪_{\varphi \theta }𝐩^T𝒜^{}𝐩^{},`$
where

$`𝒜^{}=\left(\begin{array}{ccc}c_0+\frac{4c_2m}{(1+m)^2}& c_0+\frac{2c_2}{1+m}& c_0\\ c_0+\frac{2c_2}{1+m}& c_0& c_0+\frac{2c_2}{1m}\\ c_0& c_0+\frac{2c_2}{1m}& c_0\frac{4c_2m}{(1m)^2}\end{array}\right).`$ (26)
It is then simply to verify that $`𝒜^{}`$ is semi-positive definite, which guarantees that $`𝒜`$ is positive definite. Thus the initial ground state, the mean field symmetry breaking state with a coherent condensate amplitude $`\sqrt{N}\varphi _j(\stackrel{}{r})`$ is stable. The associated quantum fluctuations of the atom numbers can be studied by the linearization approximation Eq. (3).

We note that $``$ is diagonalized by an orthogonal matrix

$`𝒰=\left(\begin{array}{ccc}1/\sqrt{3}& 1/\sqrt{3}& 1/\sqrt{3}\\ 1/\sqrt{2}& 0& 1/\sqrt{2}\\ 1/\sqrt{6}& 2/\sqrt{6}& 1/\sqrt{6}\end{array}\right),`$ (30)
which motivates the introduction of new number and phase fluctuation operators according to $`𝐩^T=𝐩^T𝒰^T=(p_N,p_M,p_Y)`$ and $`𝐪^T=𝐪^T𝒰^T=(q_N,q_M,q_Y)`$. Not surprisingly, this transformation recovers the underline symmetry of a spinor-1 condensate as discovered recently in Ref. ozgur . It turns out that $`p_N`$, $`p_M`$, and $`p_Y`$ represent respectively the fluctuations of the total number of atoms, the magnetization, and the hypercharge. $`𝒜`$ is not simultaneously diagonalized by the $`𝒰`$, yet the zero mode Hamiltonian takes a much simpler form

$`{\displaystyle \frac{H_{\mathrm{zero}}}{N}}`$ $`=`$ $`ap_N^2+bp_M^2+cp_Y^2`$ (31)
$`+\alpha p_Np_Y+\beta p_Mp_Y+\gamma p_Np_M`$
$`+\eta q_Y^2,`$
where all the coefficients are given in Appendix C. By replacing $`p_j`$ with $`i\mathrm{}/q_j`$, we obtain an eigenvalue equation for the distribution of the various fluctuations in the condensate ground state

$`{\displaystyle \frac{E}{N\mathrm{}^2}}\phi `$ $`=`$ $`a{\displaystyle \frac{^2\phi }{q_N^2}}b{\displaystyle \frac{^2\phi }{q_M^2}}c{\displaystyle \frac{^2\phi }{q_Y^2}}`$ (32)
$`\alpha {\displaystyle \frac{^2\phi }{q_Nq_Y}}\beta {\displaystyle \frac{^2\phi }{q_Mq_Y}}\gamma {\displaystyle \frac{^2\phi }{q_Nq_M}}`$
$`+{\displaystyle \frac{\eta }{\mathrm{}^2}}q_Y^2\phi .`$
We look for eigenfunctions of the form

$`\chi (q_Y)e^{i(k_Nq_N+k_Mq_M)},`$ (33)
which upon substituting into (32) yields

$`c{\displaystyle \frac{d^2\chi }{dq_Y^2}}+i\kappa {\displaystyle \frac{d\chi }{dq_Y}}+\left(E^{}{\displaystyle \frac{\eta }{\mathrm{}^2}}q_Y^2\right)\chi =0,`$ (34)
with

$`\kappa `$ $`=`$ $`\alpha k_N+\beta k_M,`$
$`E^{}`$ $`=`$ $`{\displaystyle \frac{E}{N\mathrm{}^2}}(ak_N^2+bk_M^2+\gamma k_Nk_M).`$ (35)
The eigen-equation (34) is essentially of the harmonic type, which can be easily solved and the resulting eigenvalues of Eq. (32) are

$`{\displaystyle \frac{E_n}{N}}=(2n+1)\mathrm{}(c\eta )^{1/2}+\mathrm{\Lambda }(k_N,k_M),`$
where $`n=0,1,2,\mathrm{}`$ and

$`\mathrm{\Lambda }(k_N,k_M)`$ $`=`$ $`\mathrm{}^2(ak_N^2+bk_M^2+\gamma k_Nk_M)`$
$`{\displaystyle \frac{\mathrm{}^2}{4c}}(\alpha k_N+\beta k_M)^2.`$
We note that $`\mathrm{\Lambda }(k_N,k_M)`$ must be greater than or equal to zero due to the positive definiteness of matrix $`𝒜`$. Therefore for the ground state of Eq. (32), we must have $`k_N=k_M=0`$, the ground state ($`n=0`$) energy and eigenfunction are respectively

$`E_0`$ $`=`$ $`N\mathrm{}(c\eta )^{1/2}{\displaystyle \frac{1}{2}}\mathrm{}\mathrm{\Omega },`$ (36)
$`\phi _0(q_Y)`$ $`=`$ $`\pi ^{1/4}(\mathrm{}\mathrm{\Delta }_{Y0})^{1/2}\mathrm{exp}\left({\displaystyle \frac{q_Y^2}{2\mathrm{}^2\mathrm{\Delta }_{Y0}^2}}\right),`$ (37)
with

$$\mathrm{\Omega }=\frac{2N}{\mathrm{}}[c_2𝒪_{\varphi \varphi }(\epsilon 4c_2𝒪_{\varphi \theta })]^{1/2},$$

and

$$\mathrm{\Delta }_{Y0}^2=\frac{1}{\mathrm{}}\sqrt{\frac{c}{\eta }}=\frac{1}{3n_0^2}\left(\frac{\epsilon 4c_2𝒪_{\varphi \theta }}{c_2𝒪_{\varphi \varphi }}\right)^{1/2}.$$

This ground state of $`H_{\mathrm{zero}}`$ can be easily understood. For a system starting with a fixed total number of atoms and magnetization, the dynamic conservation of the total number of atoms and total magnetization as dictated by the Hamiltonian (1) requires $`p_N^2=p_M^2=0`$, which leads to $`q_N^2=q_M^2=\mathrm{}`$, i.e., the ground state of $`H_{\mathrm{zero}}`$ has completely diffused phases for $`q_N`$ and $`q_M`$; The inherent ferromagnetic interaction, nevertheless, prepares a correlated ground state such that both $`p_Y`$ and $`q_Y`$ takes a Gaussian form distribution with $`p_Y^2=1/2\mathrm{\Delta }_{Y0}^2`$ and $`q_Y^2=\mathrm{}^2\mathrm{\Delta }_{Y0}^2/2`$, which is in fact the minimal uncertainty coherent state that is consistent with the symmetry breaking mean field assumptions of the ground state with three separate macroscopic wave functions for each of the spin components. However, it is experimentally difficult to produce a condensate having total fixed number and phase fluctuations as specified by such a state. We therefore turn to study the number and phase fluctuations of a general state. We start with the dynamic equations for $`p_j`$ and $`q_j`$,

$`\dot{p}_N`$ $`=`$ $`0,`$
$`\dot{p}_M`$ $`=`$ $`0,`$
$`\dot{p}_Y`$ $`=`$ $`2N\eta q_Y,`$
$`\dot{q}_N`$ $`=`$ $`N\left(2ap_N+\gamma p_M+\alpha p_Y\right),`$
$`\dot{q}_M`$ $`=`$ $`N\left(\gamma p_N+2bp_M+\beta p_Y\right),`$
$`\dot{q}_Y`$ $`=`$ $`N\left(\alpha p_N+\beta p_M+2cp_Y\right).`$ (38)
By noting that both $`p_Np_N(0)`$ and $`p_Mp_M(0)`$ are constants of motion, we can express the solutions as

$`𝐱(t)=𝒯(t)𝐱(0),`$ (39)
where $`𝐱^T(t)=[p_N(t),p_M(t),p_Y(t),q_N(t),q_M(t),q_Y(t)]`$, $`𝐱^T(0)=[p_N(0),p_M(0),p_Y(0),q_N(0),q_M(0),q_Y(0)]`$ is its initial value, and the evolution matrix $`𝒯(t)`$ takes following form

$`𝒯(t)=\left(\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ \alpha ^{}(1\mathrm{cos}\mathrm{\Omega }t)& \beta ^{}(1\mathrm{cos}\mathrm{\Omega }t)& \mathrm{cos}\mathrm{\Omega }t& 0& 0& \mathrm{sin}\mathrm{\Omega }t/\mathrm{}\mathrm{\Delta }_{Y0}^2\\ \zeta _aNt+\mathrm{}\mathrm{\Delta }_{Y0}^2\alpha ^2\mathrm{sin}\mathrm{\Omega }t& \zeta _cNt+\mathrm{}\mathrm{\Delta }_{Y0}^2\alpha ^{}\beta ^{}\mathrm{sin}\mathrm{\Omega }t& \mathrm{}\mathrm{\Delta }_{Y0}^2\alpha ^{}\mathrm{sin}\mathrm{\Omega }t& 1& 0& \alpha ^{}(1\mathrm{cos}\mathrm{\Omega }t)\\ \zeta _cNt+\mathrm{}\mathrm{\Delta }_{Y0}^2\alpha ^{}\beta ^{}\mathrm{sin}\mathrm{\Omega }t& \zeta _bNt+\mathrm{}\mathrm{\Delta }_{Y0}^2\beta ^2\mathrm{sin}\mathrm{\Omega }t& \mathrm{}\mathrm{\Delta }_{Y0}^2\beta ^{}\mathrm{sin}\mathrm{\Omega }t& 0& 1& \beta ^{}(1\mathrm{cos}\mathrm{\Omega }t)\\ \mathrm{}\mathrm{\Delta }_{Y0}^2\alpha ^{}\mathrm{sin}\mathrm{\Omega }t& \mathrm{}\mathrm{\Delta }_{Y0}^2\beta ^{}\mathrm{sin}\mathrm{\Omega }t& \mathrm{}\mathrm{\Delta }_{Y0}^2\mathrm{sin}\mathrm{\Omega }t& 0& 0& \mathrm{cos}\mathrm{\Omega }t\end{array}\right).`$ (46)
The coefficients $`\alpha ^{}`$, $`\beta ^{}`$, $`\zeta _a`$, $`\zeta _b`$, and $`\zeta _c`$ are given in Appendix D. The explicit solutions of the time-dependent number and phase fluctuations are given in Appendix E. We note that the general solutions has a simple structure: in addition to oscillating terms of the form $`\mathrm{cos}\mathrm{\Omega }t`$ and $`\mathrm{sin}\mathrm{\Omega }t`$, the phase fluctuations of $`q_N`$ and $`q_M`$ also contains terms proportional to $`Nt`$. This indicates that our linearization approximation Eq. (3) is valid only for a finite duration. The time-dependent covariance can also be obtained to be

$`x_i(t)x_j(t)={\displaystyle \underset{kl}{}}𝒯_{ik}(t)𝒯_{jl}(t)x_k(0)x_l(0).`$ (47)
Finally, let’s consider the diffusion of the direction of the macroscopic spin $`\stackrel{}{f}(\stackrel{}{r})_{ij}\psi _j^{}(\stackrel{}{r})\stackrel{}{𝐅}_{ij}\psi _j(\stackrel{}{r})`$. As we have shown in Ref. su , independent of the spatial coordinates, the spins of Bose condensed atoms are parallel for ferromagnetic interactions, i.e., it essentially acts as a macroscopic magnetic dipole pointing along the same direction,

$`\stackrel{}{f}_0(\stackrel{}{r})`$ $`=`$ $`N{\displaystyle \underset{ij}{}}\varphi _i^{}(\stackrel{}{r})\stackrel{}{𝐅}_{ij}\varphi _j(\stackrel{}{r})`$ (51)
$`=`$ $`N\varphi ^2(\stackrel{}{r})\left(\begin{array}{c}\sqrt{1m^2}\mathrm{cos}\mathrm{\Theta }\\ \sqrt{1m^2}\mathrm{sin}\mathrm{\Theta }\\ m\end{array}\right),`$
where $`\mathrm{\Theta }=\alpha _+\alpha _0=\alpha _0\alpha _{}`$. Our initial choice of three separate macroscopic condensate wave functions (with given phases) for each of the spin components fixes the initial direction of the condensate spin direction. As the phase spreading dynamics attempts to restore the $`U(1)`$ symmetries of each of the phases, our initial choice of phases will become irrelevant. Under the linear approximation (3) the fluctuation becomes

$`\delta \stackrel{}{f}(t)={\displaystyle 𝑑\stackrel{}{r}\left[\stackrel{}{f}(\stackrel{}{r})\stackrel{}{f}_0(\stackrel{}{r})\right]},`$ (52)
or explicitly

$`\delta f_x`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}{\displaystyle }d\stackrel{}{r}(\varphi _0^{}\delta \psi _++\varphi _+\delta \psi _0^{}+\varphi _{}^{}\delta \psi _0+\varphi _0\delta \psi _{}^{}`$
$`+h.c.),`$
$`\delta f_y`$ $`=`$ $`{\displaystyle \frac{i}{\sqrt{2}}}{\displaystyle }d\stackrel{}{r}(\varphi _0^{}\delta \psi _++\varphi _+\delta \psi _0^{}+\varphi _{}^{}\delta \psi _0+\varphi _0\delta \psi _{}^{}`$
$`h.c.),`$
$`\delta f_z`$ $`=`$ $`{\displaystyle }d\stackrel{}{r}(\varphi _+^{}\delta \psi _+\varphi _{}^{}\delta \psi _{}+h.c.).`$ (53)
Although $`\delta \stackrel{}{f}(t)0`$, the dynamics of the zero mode causes the variance of $`\delta \stackrel{}{f}(t)`$ to be nonzero. It can be reexpressed in terms of the (time-dependent) number and phase fluctuation operators under the RWA as

$`\delta f_x(t)`$ $`=`$ $`{\displaystyle \frac{\mathrm{cos}\mathrm{\Theta }}{\sqrt{1m^2}}}\left[\sqrt{3}p_N(t)\sqrt{2}mp_M(t)\right]`$
$`+{\displaystyle \frac{\sqrt{1m^2}\mathrm{sin}\mathrm{\Theta }}{2\mathrm{}}}\left[\sqrt{2}q_M(t)+\sqrt{6}mq_Y(t)\right],`$
$`\delta f_y(t)`$ $`=`$ $`{\displaystyle \frac{\mathrm{sin}\mathrm{\Theta }}{\sqrt{1m^2}}}\left[\sqrt{3}p_N(t)\sqrt{2}mp_M(t)\right]`$
$`+{\displaystyle \frac{\sqrt{1m^2}\mathrm{cos}\mathrm{\Theta }}{2\mathrm{}}}\left[\sqrt{2}q_M(t)+\sqrt{6}mq_Y(t)\right],`$
$`\delta f_z(t)`$ $`=`$ $`\sqrt{2}p_M(t),`$ (54)
where the time-dependent terms are

$`\sqrt{2}q_M(t)+\sqrt{6}mq_Y(t)=\sqrt{2}q_M(0)+\sqrt{6}mq_Y(0)`$
$`+{\displaystyle \frac{4\epsilon Nt}{1m^2}}\left[\sqrt{3}mp_N(0)+\sqrt{2}p_M(0)\right].`$ (55)
We see that $`\delta f_z`$ is a constant of motion due to the conservation of total magnetization. On the other hand, both $`\delta f_x`$ and $`\delta f_y`$ diffuse with time. To calculate the variance of $`\delta \stackrel{}{f}(t)`$, we simply take $`\mathrm{\Theta }=0`$, which yields

$`\delta f_x(t)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{1m^2}}}\left[\sqrt{3}p_N(0)\sqrt{2}mp_M(0)\right],`$
$`\delta f_y(t)`$ $`=`$ $`{\displaystyle \frac{\sqrt{1m^2}}{2\mathrm{}}}\left[\sqrt{2}q_M(t)+\sqrt{6}mq_Y(t)\right],`$
$`\delta f_z(t)`$ $`=`$ $`\sqrt{2}p_M(0),`$ (56)
i.e., with the direction of the initial macroscopic spin constrained in the x-z plane. The results of the variances for condensate spin direction fluctuations are listed in Appendix F.

We can also project the fluctuation of the spin \[Eq. (54)\] onto the spherical coordinates $`(\widehat{r},\widehat{\theta },\widehat{\varphi })`$ as illustrated in Fig. 5, and noting that

$`\widehat{r}`$ $`=`$ $`\widehat{x}\mathrm{sin}\theta \mathrm{cos}\varphi +\widehat{y}\mathrm{sin}\theta \mathrm{sin}\varphi +\widehat{z}\mathrm{cos}\theta ,`$
$`\widehat{\theta }`$ $`=`$ $`\widehat{x}\mathrm{cos}\theta \mathrm{cos}\varphi +\widehat{y}\mathrm{cos}\theta \mathrm{sin}\varphi \widehat{z}\mathrm{sin}\theta ,`$
$`\widehat{\varphi }`$ $`=`$ $`\widehat{x}\mathrm{sin}\varphi +\widehat{y}\mathrm{cos}\varphi ,`$ (57)
and

$`\mathrm{sin}\theta =\sqrt{1m^2},`$ $`\mathrm{cos}\theta =m,`$
$`\mathrm{sin}\varphi =\mathrm{sin}\mathrm{\Theta },`$ $`\mathrm{cos}\varphi =\mathrm{cos}\mathrm{\Theta },`$
we find

$`\delta f_r(t)`$ $`=`$ $`\sqrt{3}p_N(0),`$
$`\delta f_\theta (t)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{1m^2}}}\left[\sqrt{3}mp_N(0)\sqrt{2}p_M(0)\right],`$
$`\delta f_\varphi (t)`$ $`=`$ $`{\displaystyle \frac{\sqrt{1m^2}}{2\mathrm{}}}\left[\sqrt{2}q_M(t)+\sqrt{6}mq_Y(t)\right].`$ (58)
We immediately see that both $`\delta f_r(t)\delta f_r(0)`$ and $`\delta f_\theta (t)\delta f_\theta (0)`$ are fixed due to the conservations of $`N`$ and $``$. Furthermore, with the use of Eq. (55), we find that $`\delta f_\varphi (t)`$ can be rewritten as

$`\delta f_\varphi (t)=\delta f_\varphi (0){\displaystyle \frac{2\epsilon Nt}{\mathrm{}}}\delta f_\theta (0),`$ (59)
where

$$\delta f_\varphi (0)=\frac{\sqrt{1m^2}}{2\mathrm{}}\left[\sqrt{2}q_M(0)+\sqrt{6}mq_Y(0)\right].$$

We then see that $`\delta f_\varphi (t)`$ grows linearly with $`t`$ and the diffusion rate is given by

$`R_d={\displaystyle \frac{2\epsilon N}{\mathrm{}}}\delta f_\theta (0),`$ (60)
which is proportional to $`\delta f_\theta (0)`$. A result easily understandable in terms of the single axis twisting of the isospin $`T_z^2`$ (of the Hamiltonian) of a spior-1 spinor condensate ozgur .

Given any initial condensate state, its variances in atom numbers for each individual spin component and their correlations and related phase fluctuations completely determine the subsequent fluctuations of $`\delta f_\varphi (t)`$. $`P_j(0)`$ is in fact nothing but the atom number fluctuation of state $`|j`$. In general we we note that

$`{\displaystyle 𝑑\stackrel{}{r}\psi _j^{}(\stackrel{}{r})\psi _j(\stackrel{}{r})}`$ $`=`$ $`N+\sqrt{N}P_j+{\displaystyle 𝑑\stackrel{}{r}\delta \psi _j^{}(\stackrel{}{r})\delta \psi _j(\stackrel{}{r})}`$ (61)
$`=`$ $`N+\sqrt{N}P_j`$
$`+{\displaystyle \frac{sn_j}{2J_j^2}}P_j^2+{\displaystyle \frac{n_j}{2\mathrm{}^2J_j^2}}Q_j^2{\displaystyle \frac{1}{2}},`$
where $`s=𝑑\stackrel{}{r}\theta ^2(\stackrel{}{r})`$ and we have used RWA in deriving the last line. One can easily find

$`\delta N_j^2`$ $`=`$ $`N_j^2N_j^2`$ (62)
$`=`$ $`\left[{\displaystyle 𝑑\stackrel{}{r}\psi _j^{}(\stackrel{}{r})\psi _j(\stackrel{}{r})}\right]^2{\displaystyle 𝑑\stackrel{}{r}\psi _j^{}(\stackrel{}{r})\psi _j(\stackrel{}{r})}^2`$
$``$ $`NP_j^2=Np_j^2,`$
to lowest order.

If we use $`\sigma _N^2`$, $`\sigma _M^2`$, and $`\sigma _Y^2`$ to denote the variances of the total atom number ($`N`$), the magnetization ($``$), and the hyper-charge ($`Y`$), we find that

$`\left[\delta f_r(t)\right]^2`$ $`=`$ $`3\sigma _N^2,`$
$`\left[\delta f_\theta (t)\right]^2`$ $`=`$ $`{\displaystyle \frac{1}{1m^2}}\left(3m^2\sigma _N^2+2\sigma _M^2\right).`$ (63)
To calculate the variance of $`\delta f_\varphi (t)`$, we assume a completely uncorrelated (between fluctuations of $`N`$, $``$, and $`Y`$) initial distribution in the form of a Gaussian wave packet

$`\phi (p_N,p_M,p_Y)=(2\pi )^{3/4}(\sigma _N\sigma _M\sigma _Y)^{1/2}`$
$`\times \mathrm{exp}\left[{\displaystyle \frac{1}{4}}\left({\displaystyle \frac{p_N^2}{\sigma _N^2}}+{\displaystyle \frac{p_M^2}{\sigma _M^2}}+{\displaystyle \frac{p_Y^2}{\sigma _Y^2}}\right)\right].`$ (64)
Since $`q_j=i\mathrm{}/p_j`$, we find for such a distribution

$`q_i(0)q_j(0)`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{4\sigma _i^2}}\delta _{ij},`$
$`p_i(0)p_j(0)`$ $`=`$ $`\sigma _i^2\delta _{ij},`$
$`q_i(0)p_i(0)+p_i(0)q_i(0)`$ $`=`$ $`0,`$
$`q_i(0)p_j(0)`$ $`=`$ $`0,\text{ for }ij,`$ (65)
and $`i,j=N,M,Y`$. Therefore

$`\left[\delta f_\varphi (t)\right]^2`$ $`=`$ $`\left[\delta f_\varphi (0)\right]^2+{\displaystyle \frac{4\epsilon ^2N^2t^2}{\mathrm{}^2}}\left[\delta f_\theta (0)\right]^2,`$ (66)
with

$`\left[\delta f_\varphi (0)\right]^2={\displaystyle \frac{1m^2}{8}}\left({\displaystyle \frac{1}{\sigma _M^2}}+3m^2{\displaystyle \frac{1}{\sigma _Y^2}}\right).`$ (67)
If on the other hand, we take a distribution similar to Eq. (64) but with no correlations between the populations in the spin components $`j=+,0,`$,

$`\phi (p_+^{},p_0^{},p_{}^{})=(2\pi )^{3/4}(\sigma _+\sigma _0\sigma _{})^{1/2}`$
$`\times \mathrm{exp}\left[{\displaystyle \frac{1}{4}}\left({\displaystyle \frac{p_+^22}{\sigma _+^2}}+{\displaystyle \frac{p_0^2}{\sigma _0^2}}+{\displaystyle \frac{p_{}^2}{\sigma _{}^2}}\right)\right],`$ (68)
with $`\sigma _j^2`$ denoting the corresponding atom number variance, we find

$`\left[\delta f_r(t)\right]^2`$ $`=`$ $`\sigma _+^2+\sigma _0^2+\sigma _{}^2,`$
$`\left[\delta f_\theta (t)\right]^2`$ $`=`$ $`{\displaystyle \frac{1}{1m^2}}[(1m)^2\sigma _+^2+m^2\sigma _0^2`$ (69)
$`+(1+m)^2\sigma _{}^2],`$
and Eq. (66) still holds but with

$`\left[\delta f_\varphi (0)\right]^2`$ $`=`$ $`{\displaystyle \frac{1m^2}{16}}[{\displaystyle \frac{(1+m)^2}{\sigma _+^2}}+4m^2{\displaystyle \frac{1}{\sigma _0^2}}`$ (70)
$`+{\displaystyle \frac{(1m)^2}{\sigma _{}^2}}].`$
In the Thomas-Fermi (TF) regime, we obtain most of the above results analytically. First we note that the ground wave function is

$`\varphi ^2(\stackrel{}{r})=\{\begin{array}{cc}(c_+N)^1[\mu V_{\mathrm{ext}}(\stackrel{}{r})],& \text{if }\mu >V_{\mathrm{ext}}(\stackrel{}{r}),\hfill \\ 0,& \text{otherwise.}\hfill \end{array},`$ (73)
where $`c_+=4\pi \mathrm{}^2a_2/M`$ and the chemical potential

$$\mu =\frac{\mathrm{}\omega _{\mathrm{ho}}}{2}\left(\frac{15Na_2}{a_{\mathrm{ho}}}\right)^{2/5}$$

with $`\omega _{\mathrm{ho}}=(\omega _x\omega _y\omega _z)^{1/3}`$ and $`a_{\mathrm{ho}}=(\mathrm{}/M\omega _{\mathrm{ho}})^{1/2}`$. Similarly, ignoring the kinetic energy term in the Eq. (18), we obtain

$`\theta (\stackrel{}{r})=\{\begin{array}{cc}[2V_0\varphi (\stackrel{}{r})]^1,& \text{if }\mu >V_{\mathrm{ext}}(\stackrel{}{r}),\\ 0,& \text{otherwise},\end{array}`$ (76)
and

$`\stackrel{~}{u}={\displaystyle \frac{c_+}{V_0}}`$ (77)
with

$$V_0=_{\mu V_{\mathrm{ext}}(\stackrel{}{r})>0}𝑑\stackrel{}{r}=\frac{4\pi a_{\mathrm{ho}}^3}{3}\left(\frac{15a_2}{a_{\mathrm{ho}}}\right)^{3/5}N^{3/5}.$$

All overlap integrals can be calculated analytically as

$`𝒪_{\varphi \theta }`$ $`=`$ $`(4V_0)^1={\displaystyle \frac{3}{16\pi a_{\mathrm{ho}}^3}}\left({\displaystyle \frac{a_{\mathrm{ho}}}{15a_2}}\right)^{3/5}N^{3/5},`$
$`𝒪_{\varphi \varphi }`$ $`=`$ $`{\displaystyle \frac{32\pi \mu ^{7/2}}{105c_+^2N^2\omega _{\mathrm{ho}}^3}}\left({\displaystyle \frac{2}{M}}\right)^{3/2}`$
$`=`$ $`{\displaystyle \frac{15^{2/5}}{14\pi a_{\mathrm{ho}}^3}}\left({\displaystyle \frac{a_{\mathrm{ho}}}{a_2}}\right)^{3/5}N^{3/5}.`$
The parameter $`\epsilon `$ relates essentially to the kinetic energy operator $`E_k^{(\theta )}=(\mathrm{}^2/2M)𝑑\stackrel{}{r}\theta ^2\theta `$ evaluated with respect to the phase function $`\theta (\stackrel{}{r})`$. In the TF approximation we find

$`N\epsilon =2E_k^{(\theta )}.`$ (78)
Using Eqs. (15) and (18) we find

$`(\mu +3c_+N\varphi ^2)\theta ={\displaystyle \frac{N\stackrel{~}{u}}{\mu }}(+c_+N\varphi ^2)\varphi ,`$ (79)
which simplifies to

$`\left({\displaystyle \frac{\mathrm{}^2^2}{2M}}+2c_+N\varphi ^2\right)\theta ={\displaystyle \frac{N\stackrel{~}{u}}{\mu }}\left({\displaystyle \frac{\mathrm{}^2^2}{2M}}+\mu \right)\varphi ,`$ (80)
in the TF limit. Multiplying this Eq. (80) from left by $`\theta `$ or $`\varphi `$ respectively, and integrating over $`\stackrel{}{r}`$, we obtain

$`E_k^{(\theta )}+2c_+N𝒪_{\varphi \theta }`$ $`=`$ $`{\displaystyle \frac{N\stackrel{~}{u}}{\mu }}\left(E_k^{(\varphi \theta )}+{\displaystyle \frac{\mu }{2}}\right),`$
$`E_k^{(\varphi \theta )}+2c_+N{\displaystyle 𝑑\stackrel{}{r}\varphi ^3\theta }`$ $`=`$ $`{\displaystyle \frac{N\stackrel{~}{u}}{\mu }}\left(E_k^{(\varphi )}+\mu \right),`$ (81)
where $`E_k^{(\varphi )}=(\mathrm{}^2/2M)𝑑\stackrel{}{r}\varphi ^2\varphi `$ is the kinetic energy per atom in the condensate, and $`E_k^{(\varphi \theta )}=(\mathrm{}^2/2M)𝑑\stackrel{}{r}\varphi ^2\theta =(\mathrm{}^2/2M)𝑑\stackrel{}{r}\theta ^2\varphi `$. Eliminating $`E_k^{(\varphi \theta )}`$ from above two equations, we then obtain

$`E_k^{(\theta )}`$ $`=`$ $`\left({\displaystyle \frac{N\stackrel{~}{u}}{\mu }}\right)^2E_k^{(\varphi )}+{\displaystyle \frac{N\stackrel{~}{u}}{2}}2c_+N𝒪_{\varphi \theta }`$ (82)
$`+{\displaystyle \frac{N\stackrel{~}{u}}{\mu }}\left(N\stackrel{~}{u}2c_+N{\displaystyle 𝑑\stackrel{}{r}\varphi ^3\theta }\right)`$
$`=`$ $`\left({\displaystyle \frac{N\stackrel{~}{u}}{\mu }}\right)^2E_k^{(\varphi )}`$
$`=`$ $`{\displaystyle \frac{4}{25}}E_k^{(\varphi )}.`$
For a spherical trap, adopting the result of Ref. dalfovo , we immediately find

$`\epsilon ={\displaystyle \frac{4\mathrm{}^2}{5MNR^2}}\mathrm{ln}\left({\displaystyle \frac{R}{1.3a_{\mathrm{ho}}}}\right),`$ (83)
with

$$R=(15a_2a_{\mathrm{ho}}^4)^{1/5}N^{1/5},$$

which seems to indicate that $`N\epsilon N^{2/5}`$, asymptotically goes to zero at large $`N`$. By all measures, we find $`\epsilon `$ to be a very small quantity as compared to $`c_2𝒪_{\varphi \theta }`$. Taking $`\epsilon `$ as zero in the TF limit, we obtain

$`n_0^2\mathrm{\Delta }_{Y0}^2`$ $`=`$ $`{\displaystyle \frac{1}{3}}\sqrt{{\displaystyle \frac{7}{10}}},`$
$`\mathrm{}\mathrm{\Omega }`$ $`=`$ $`{\displaystyle \frac{\sqrt{3/14}}{15^{1/10}}}\left({\displaystyle \frac{c_2}{\pi a_{\mathrm{ho}}^3}}\right)\left({\displaystyle \frac{a_{\mathrm{ho}}}{a_2}}\right)^{3/5}N^{2/5}.`$ (84)
To shed light on the above scaling results based on the TF approximation, we have performed extensive numerical calculations. As an example, we consider a spinor condensate of <sup>87</sup>Rb atoms in a spherically symmetric trap with a harmonic frequency $`\omega =(2\pi )100`$ (Hz) and consider the regime of $`N10^710^{15}`$. It turns out that our results indeed confirm the scaling relationship

$$𝒪_{\varphi \varphi },𝒪_{\varphi \theta },\stackrel{~}{u},\mathrm{\Omega }/NN^{3/5},$$

and

$`{\displaystyle \frac{N\stackrel{~}{u}}{\mu }}`$ $``$ $`{\displaystyle \frac{2}{5}},`$
$`n_0^2\mathrm{\Delta }_{Y0}^2`$ $``$ $`{\displaystyle \frac{1}{3}}\sqrt{{\displaystyle \frac{7}{10}}},`$ (85)
in the large $`N`$ limit. We also find that

$$E_k^{(\varphi )}N^{0.3549},$$

in rough agreement with the result of Ref. dalfovo ; On the other hand, we find

$$N\epsilon N^{0.1348},E_k^{(\theta )}N^{0.1314},$$

and

$$N\epsilon 4E_k^{(\theta )},$$

which do not agree with the results above as obtained with the TF approximation. A careful analysis as confirmed by our numerical results reveal that although $`𝑑\stackrel{}{r}(V_{\mathrm{ext}}\mu +c_+N\varphi ^2)\varphi ^2`$ is small, $`𝑑\stackrel{}{r}(V_{\mathrm{ext}}\mu +c_+N\varphi ^2)\theta ^2`$ is not, and approximately we find

$`E_k^{(\theta )}{\displaystyle \frac{1}{2}}{\displaystyle 𝑑\stackrel{}{r}(V_{\mathrm{ext}}\mu +c_+N\varphi ^2)\theta ^2}.`$ (86)
This leads to the opposite $`N`$-dependence for $`E_k^{(\theta )}`$ and $`E_k^{(\varphi )}`$; while $`E_k^{(\theta )}`$ increases with increasing values of $`N`$, $`E_k^{(\varphi )}`$ decreases. Thus the identity Eq. (82) as obtained in the TF limit is invalid.

## IV A nonzero magnetic field

Inside a nonzero homogeneous magnetic field ($`\stackrel{}{B}`$), the atomic quantization axis becomes fixed, along the direction of $`\stackrel{}{B}`$, or more conveniently denoted as the z-axis. Mathematically, it turns out this is equivalent to the the requirement of the conservation of magnetization $`=N_+N_{}`$. This reduces the SO(3) rotational symmetry to its subgroup SO(2), which is still a continuous symmetry. Thus the quantum dynamics of the phase remains formally the same as discussed before.

In the linear Zeeman regime, the original Hamiltonian (1) is augmented to

$`H^{}=H\mathrm{}\omega _L𝐅_z.`$ (87)
The Larmor precessing frequency is $`\omega _L=B\mu _B/\mathrm{}`$ with $`\mu _B`$ the magnetic dipole moment for state $`|F=1,M_F=1`$. In this case, the only change is the phases of the spatial wave functions and the phase functions if of the following time-dependent form

$`\alpha _j=\alpha _{j0}+j\omega _Lt,`$ (88)
for $`j=+,0,`$ and $``$. $`\alpha _{+0}\alpha _{00}=\alpha _{00}\alpha _0=0`$ still holds true for a ferromagnetic spinor-1 condensate. A direct consequence of the magnetic field is thus a time-dependent relative phase

$`\mathrm{\Theta }(t)`$ $`=`$ $`\alpha _{+0}\alpha _{00}+\omega _Lt`$ (89)
$`=`$ $`\alpha _{00}\alpha _0+\omega _Lt.`$
When applied to the macroscopic condensate spin, we find the spin gains a precession around the z-axis.

Now the picture of the dynamics of the macroscopic condensate spin is clear: besides the precession induced by the external B-field, it also diffuses to restore the U(1) and SO(2) symmetries of the Hamiltonian.

The effect of quadratic Zeeman shift can also be simply addressed. Apart from an overall shift, it introduces a (positive) B-field dependent level shift $`\mathrm{}\delta _BF_z^2`$ to states $`|\pm `$ with respect to state $`|0`$. This differential shift causes the precessing of the condensate spin to be twisted, i.e. the precessing of the $`|+`$ ($`|`$) component being slower (faster) by $`\delta _B`$.

## V The case of anti-ferromagnetic interactions

Our discussion of the phase diffusion in the previous section is based on the proof that the single spatial mode approximation is exact for ferromagnetic interactions. One might therefore conclude that a similar study for anti-ferromagnetic interactions would be difficult as, unlike for ferromagnetic interactions, it’s difficult to predict the ground state condensate structure. In reality, however, the quantum phase problem for anti-ferromagnetic interactions is much simpler, provided the mean field description for the ground state still applies.

First we note that when $`0`$, the population in the $`0`$-th component is zero in the ground state pu ; su . Therefore, the zero mode Hamiltonian (20) reduces to

$`{\displaystyle \frac{H_{\mathrm{zero}}}{N}}`$ $`=`$ $`{\displaystyle \frac{P_+^2}{J_+^2}}\left({\displaystyle \frac{J_+}{2}}\stackrel{~}{u}+2c_{}I_+\right)`$ (90)
$`{\displaystyle \frac{P_{}^2}{J_{}^2}}\left({\displaystyle \frac{J_{}}{2}}\stackrel{~}{u}+2c_{}I_+\right)`$
$`+{\displaystyle \frac{4c_{}P_+P_{}}{J_+J_{}}}I_+,`$
which is exactly the same zero mode Hamiltonian for a binary condensate with no coupling between its two components. This has already been studied before villain97 . We note that results obtained in Section II remain valid.

Now, let’s consider the $`=0`$ case. As was pointed out in Ref. su , the single mode approximation again applies in this case, and we can assume

$`\varphi _j(\stackrel{}{r})=\sqrt{n_j}\varphi (\stackrel{}{r})e^{i\alpha _j},`$ (91)
where $`n_+=n_{}=(1n_0)/2`$ and total energy of the system is independent of the value of $`n_0`$. A similar relation exists among the phases of the three components: $`\alpha _++\alpha _{}2\alpha _0=\pi `$. The spatial profile $`\varphi (\stackrel{}{r})`$ satisfies following equation

$`\left[+c_0N\varphi ^2(\stackrel{}{r})\right]\varphi (\stackrel{}{r})=\mu \varphi (\stackrel{}{r}),`$ (92)
with a coefficient $`c_0`$ for the nonlinear term. The phase function then becomes

$`\theta _j(\stackrel{}{r})=\sqrt{n_j}\theta (\stackrel{}{r})e^{i\alpha _j},`$ (93)
where $`\theta (\stackrel{}{r})`$ is the solution of the following equation

$`\left[\mu +3c_0N\varphi ^2(\stackrel{}{r})\right]\theta (\stackrel{}{r})=N\stackrel{~}{u}\theta (\stackrel{}{r}).`$ (94)
Following the results of Sect. III, we define

$`\epsilon ={\displaystyle \frac{2}{N}}{\displaystyle 𝑑\stackrel{}{r}\theta (\stackrel{}{r})\left[\mu +c_0N\varphi ^2(\stackrel{}{r})\right]\theta (\stackrel{}{r})},`$ (95)
and the zero mode Hamiltonian becomes

$`{\displaystyle \frac{H_{\mathrm{zero}}}{N}}`$ $`=`$ $`\left[{\displaystyle \frac{\epsilon }{(1n_0)𝒪_{\varphi \theta }}}+2c_0+{\displaystyle \frac{2c_2}{1n_0}}\right]p_+^2𝒪_{\varphi \theta }`$ (96)
$`+\left({\displaystyle \frac{\epsilon }{2n_0𝒪_{\varphi \theta }}}+2c_0\right)p_0^2𝒪_{\varphi \theta }`$
$`+\left[{\displaystyle \frac{\epsilon }{(1n_0)𝒪_{\varphi \theta }}}+2c_0+{\displaystyle \frac{2c_2}{1n_0}}\right]p_{}^2𝒪_{\varphi \theta }`$
$`+4c_0p_+^{}p_0^{}𝒪_{\varphi \theta }+4c_0p_{}^{}p_0^{}𝒪_{\varphi \theta }`$
$`+4\left(c_0+c_2{\displaystyle \frac{2n_01}{1n_0}}\right)p_+^{}p_{}^{}𝒪_{\varphi \theta }`$
$`+{\displaystyle \frac{c_2n_0(1n_0)𝒪_{\varphi \varphi }}{2\mathrm{}^2}}\left(q_+^{}2q_0^{}+q_{}^{}\right)^2`$
$`=`$ $`𝐩^T𝒜𝐩^{}+𝐪^T𝐪^{},`$
the positive definiteness of matrices $`𝒜`$ and $``$ can again be verified. By applying the same transformation (30), we obtain

$`{\displaystyle \frac{H_{\mathrm{zero}}}{N}}`$ $`=`$ $`ap_N^2+bp_M^2+cp_Y^2+\alpha p_Np_Y+\eta q_Y^2,`$ (97)
where all coefficients are now given in Appendix G. We see that fluctuations in $`N`$ are completely decoupled from that of $`M`$ and $`Y`$ in this case, a result that also happens for ferromagnetic interactions when $`=0`$. The ground state of (96) is similar to Eq. (20) except now that

$`\mathrm{\Omega }`$ $`=`$ $`2N\sqrt{c\eta }={\displaystyle \frac{2N}{\mathrm{}}}[c_2𝒪_{\varphi \varphi }(\epsilon +4c_2n_0^2𝒪_{\varphi \theta })]^{1/2},`$
$`\mathrm{\Delta }_{Y0}^2`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{}}}\sqrt{{\displaystyle \frac{c}{\eta }}}={\displaystyle \frac{1}{3}}\left[{\displaystyle \frac{\epsilon +4c_2n_0^2𝒪_{\varphi \theta }}{c_2n_0^2(1n_0)^2𝒪_{\varphi \varphi }}}\right]^{1/2}.`$ (98)
The dynamic equations for $`p_i`$ and $`q_i`$ then become

$`\dot{p}_N`$ $`=`$ $`0,`$
$`\dot{p}_M`$ $`=`$ $`0,`$
$`\dot{p}_Y`$ $`=`$ $`2N\eta q_Y,`$
$`\dot{q}_N`$ $`=`$ $`N\left(2ap_N+\alpha p_Y\right),`$
$`\dot{q}_M`$ $`=`$ $`2Nbp_M,`$
$`\dot{q}_Y`$ $`=`$ $`N\left(\alpha p_N+2cp_Y\right).`$ (99)
One can easily find the evolution matrix to be

$`𝒯(t)=\left(\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ \alpha ^{}(1\mathrm{cos}\mathrm{\Omega }t)& 0& \mathrm{cos}\mathrm{\Omega }t& 0& 0& \mathrm{sin}\mathrm{\Omega }t/\mathrm{}\mathrm{\Delta }_{Y0}^2\\ \zeta _aNt+\mathrm{}\mathrm{\Delta }_{Y0}^2\alpha ^2\mathrm{sin}\mathrm{\Omega }t& 0& \mathrm{}\mathrm{\Delta }_{Y0}^2\alpha ^{}\mathrm{sin}\mathrm{\Omega }t& 1& 0& \alpha ^{}(1\mathrm{cos}\mathrm{\Omega }t)\\ 0& 2bNt& 0& 0& 1& 0\\ \mathrm{}\mathrm{\Delta }_{Y0}^2\alpha ^{}\mathrm{sin}\mathrm{\Omega }t& 0& \mathrm{}\mathrm{\Delta }_{Y0}^2\mathrm{sin}\mathrm{\Omega }t& 0& 0& \mathrm{cos}\mathrm{\Omega }t\end{array}\right),`$ (106)
where

$`\alpha ^{}`$ $`=`$ $`{\displaystyle \frac{\alpha }{2c}}={\displaystyle \frac{(3n_01)\epsilon +8c_2n_0^2𝒪_{\varphi \theta }}{\sqrt{2}(\epsilon +4c_2n_0^2𝒪_{\varphi \theta })}},`$
$`\zeta _a`$ $`=`$ $`{\displaystyle \frac{4ac\alpha ^2}{2c}}=12c_0𝒪_{\varphi \theta }+{\displaystyle \frac{3\epsilon (\epsilon +4c_2n_0𝒪_{\varphi \theta })}{\epsilon +4c_2n_0^2𝒪_{\varphi \theta }}}.`$
As for $`\delta \stackrel{}{f}(t)`$ we find

$`\delta f_x`$ $`=`$ $`\mathrm{cos}\mathrm{\Theta }\sqrt{{\displaystyle \frac{2n_0}{1n_0}}}p_M+\mathrm{sin}\mathrm{\Theta }{\displaystyle \frac{\sqrt{6n_0(1n_0)}}{\mathrm{}}}q_Y,`$
$`\delta f_y`$ $`=`$ $`\mathrm{sin}\mathrm{\Theta }\sqrt{{\displaystyle \frac{2n_0}{1n_0}}}p_M+\mathrm{cos}\mathrm{\Theta }{\displaystyle \frac{\sqrt{6n_0(1n_0)}}{\mathrm{}}}q_Y,`$
$`\delta f_z`$ $`=`$ $`\sqrt{2}p_M.`$ (107)
## VI conclusions

In conclusion, we have investigated in detail quantum phase diffusions of a spinor-1 condensate. When the elastic atomic interaction is of ferromagnetic type, the structure of the ground state is greatly simplified: the spatial mode functions of the different spin components are exactly the same and accordingly we also proved that the phase functions are exactly the same. This simplification allows us to construct the zero mode Hamiltonian which describes the number and phase fluctuations of a condensate while maintaining the conservations of the total number of atoms and the total magnetization. We have provided analytical results for the number and phase fluctuations due to both the quantum phase diffusion dynamics and the initial distribution of atom number and phase fluctuations. The structure of these fluctuations is rather simple: along with the oscillations due to $`\mathrm{sin}\mathrm{\Omega }t`$ and $`\mathrm{cos}\mathrm{\Omega }t`$ terms, quantum phase diffusion terms proportional to $`Nt`$ also exist. This suggests that the mean-field approach to spinor-1 condensates is only valid for a finite duration. Based on these results, we have studied the diffusion of the direction of a macroscopic condensate spin. Our investigation sheds important light on the studies of beyond mean field theory quantum correlations among Bose condensed atoms.

## VII Acknowledgements

We thank Yueheng Lan and Dr. X. X. Yi for helpful discussions. This work is supported by NSF and by a grant from the National Security Agency (NSA), Advanced Research and Development Activity (ARDA), and the Defense Advanced Research Projects Agency (DARPA) under Army Research Office (ARO) Contract No. DAAD19-01-1-0667.

## Appendix A overlap integrals

In this appendix, we list the overlap integrals as involved in the dynamic Eq. (12). Using the phase conventions as introduced in Eqs. (14) and (19), we find explicitly

$`J_j`$ $`=`$ $`n_j,`$
$`I_{jk}`$ $`=`$ $`n_jn_k𝒪_{\varphi \theta },`$
$`(\varphi _0\theta _0|\varphi _{}\theta _+)`$ $`=`$ $`(\varphi _0\theta _0|\varphi _+\theta _{})`$
$`=`$ $`(\varphi _0^2|\theta _+\theta _{})`$
$`=`$ $`n_0^2𝒪_{\varphi \theta }/2`$
$`=`$ $`2n_+n_{}𝒪_{\varphi \theta },`$
$`(\varphi _0^2|\varphi _+\varphi _{})`$ $`=`$ $`n_0^2𝒪_{\varphi \varphi }/2.`$ (108)
## Appendix B The zero mode Hamiltonian

The zero mode Hamiltonian is listed below.

$`{\displaystyle \frac{H_{\mathrm{zero}}}{N}}`$ $`=`$ $`{\displaystyle \frac{P_+^2}{J_+^2}}\left[{\displaystyle \frac{J_+}{2}}\stackrel{~}{u}+2c_{}I_++2c_+I_{+0}+2c_2(\varphi _0\theta _0|\varphi _{}\theta _+)^{}+c_2(\varphi _0^2|\theta _+\theta _{})^{}\right]`$ (109)
$`{\displaystyle \frac{2P_0^2}{J_0^2}}\left[{\displaystyle \frac{J_0}{4}}\stackrel{~}{u}+c_+I_{+0}+c_+I_0+c_2(\varphi _{}\theta _+|\varphi _0\theta _0)^{}+c_2(\varphi _+\theta _{}|\varphi _0\theta _0)^{}\right]`$
$`{\displaystyle \frac{P_{}^2}{J_{}^2}}\left[{\displaystyle \frac{J_{}}{2}}\stackrel{~}{u}+2c_{}I_++2c_+I_0+2c_2(\varphi _0\theta _0|\varphi _+\theta _{})^{}+c_2(\varphi _0^2|\theta _+\theta _{})^{}\right]`$
$`+{\displaystyle \frac{4P_+P_0}{J_+J_0}}\left[c_+I_{+0}+c_2(\varphi _0\theta _0|\varphi _{}\theta _+)^{}\right]+{\displaystyle \frac{4P_{}P_0}{J_{}J_0}}\left[c_+I_0+c_2(\varphi _0\theta _0|\varphi _+\theta _{})^{}\right]`$
$`+{\displaystyle \frac{2P_+P_{}}{J_+J_{}}}\left[2c_{}I_++c_2(\varphi _0^2|\theta _+\theta _{})^{}\right]{\displaystyle \frac{c_2}{\mathrm{}^2}}(\varphi _0^2|\varphi _+\varphi _{})^{}\left({\displaystyle \frac{Q_+}{J_+}}{\displaystyle \frac{2Q_0}{J_0}}+{\displaystyle \frac{Q_{}}{J_{}}}\right)^2.`$
## Appendix C Coefficients in the zero mode Hamiltonian for ferromagnetic interactions

The coefficients in the zero mode Hamiltonian (31) are

$`a`$ $`=`$ $`{\displaystyle \frac{1}{4n_0^2}}\mathbf{[}{\displaystyle \frac{\epsilon }{3}}(5+3m^2)`$
$`+2𝒪_{\varphi \theta }(3c_0(1m^2)^2+{\displaystyle \frac{8}{3}}c_2(13m^2))\mathbf{]},`$
$`b`$ $`=`$ $`{\displaystyle \frac{1}{2n_0^2}}\left[\epsilon (1+m^2)8c_2𝒪_{\varphi \theta }m^2\right]`$
$`c`$ $`=`$ $`{\displaystyle \frac{1}{3n_0^2}}(\epsilon 4c_2𝒪_{\varphi \theta }),`$
$`\alpha `$ $`=`$ $`{\displaystyle \frac{\sqrt{2}}{6n_0^2}}(1+3m^2)(\epsilon 4c_2𝒪_{\varphi \theta }),`$
$`\beta `$ $`=`$ $`{\displaystyle \frac{2m}{\sqrt{3}n_0^2}}(\epsilon 4c_2𝒪_{\varphi \theta }),`$
$`\gamma `$ $`=`$ $`{\displaystyle \frac{2m}{n_0^2}}\sqrt{{\displaystyle \frac{2}{3}}}\left[\epsilon c_2𝒪_{\varphi \theta }(1+3m^2)\right],`$
$`\eta `$ $`=`$ $`3c_2n_0^2𝒪_{\varphi \varphi }/\mathrm{}^2.`$ (110)
## Appendix D Parameters in Evolution matrix

The list of parameters in the evolution matrix $`𝒯(t)`$ (46) are

$`\alpha ^{}`$ $`=`$ $`{\displaystyle \frac{\alpha }{2c}}={\displaystyle \frac{1+3m^2}{2\sqrt{2}}},`$
$`\beta ^{}`$ $`=`$ $`{\displaystyle \frac{\beta }{2c}}=\sqrt{3}m,`$
$`\zeta _a`$ $`=`$ $`{\displaystyle \frac{4ac\alpha ^2}{2c}}={\displaystyle \frac{3(1+m^2)\epsilon }{1m^2}}+12(c_0+c_2)𝒪_{\varphi \theta },`$
$`\zeta _b`$ $`=`$ $`{\displaystyle \frac{4bc\beta ^2}{2c}}={\displaystyle \frac{4\epsilon }{1m^2}},`$
$`\zeta _c`$ $`=`$ $`{\displaystyle \frac{2c\gamma \alpha \beta }{2c}}={\displaystyle \frac{2\sqrt{6}m\epsilon }{1m^2}}.`$ (111)
## Appendix E Explicit form of number and phase fluctuation

Here we present the explicit form of the solutions for the number and phase fluctuations.

$`p_N(t)`$ $`=`$ $`p_N(0),`$
$`p_M(t)`$ $`=`$ $`p_M(0),`$
$`p_Y(t)`$ $`=`$ $`{\displaystyle \frac{1+3m^2}{2\sqrt{2}}}(1\mathrm{cos}\mathrm{\Omega }t)p_N(0)+\sqrt{3}m(1\mathrm{cos}\mathrm{\Omega }t)p_M(0)+p_Y(0)\mathrm{cos}\mathrm{\Omega }t{\displaystyle \frac{q_Y(0)}{\mathrm{}\mathrm{\Delta }_{Y0}^2}}\mathrm{sin}\mathrm{\Omega }t,`$
$`q_N(t)`$ $`=`$ $`\left[Nt\left({\displaystyle \frac{3(1+m^2)\epsilon }{1m^2}}+12(c_0+c_2)𝒪_{\varphi \theta }\right)+{\displaystyle \frac{(1+3m^2)^2}{8}}\mathrm{}\mathrm{\Delta }_{Y0}^2\mathrm{sin}\mathrm{\Omega }t\right]p_N(0)`$
$`\left[Nt{\displaystyle \frac{2\sqrt{6}m\epsilon }{1m^2}}+{\displaystyle \frac{\sqrt{3}m(1+3m^2)}{2\sqrt{2}}}\mathrm{}\mathrm{\Delta }_{Y0}^2\mathrm{sin}\mathrm{\Omega }t\right]p_M(0)`$
$`+{\displaystyle \frac{1+3m^2}{2\sqrt{2}}}\mathrm{}^2\mathrm{\Delta }_{Y0}^2p_Y(0)\mathrm{sin}\mathrm{\Omega }t+q_N(0){\displaystyle \frac{1+3m^2}{2\sqrt{2}}}(1\mathrm{cos}\mathrm{\Omega }t)q_Y(0),`$
$`q_M(t)`$ $`=`$ $`\left[Nt{\displaystyle \frac{2\sqrt{6}m\epsilon }{1m^2}}+\mathrm{}\mathrm{\Delta }_{Y0}^2{\displaystyle \frac{\sqrt{3}m(1+3m^2)}{2\sqrt{2}}}\mathrm{sin}\mathrm{\Omega }t\right]p_N(0)`$
$`+\left(Nt{\displaystyle \frac{4\epsilon }{1m^2}}+3m^2\mathrm{}\mathrm{\Delta }_{Y0}^2\mathrm{sin}\mathrm{\Omega }t\right)p_M(0)\sqrt{3}m\mathrm{}\mathrm{\Delta }_{Y0}^2p_Y(0)\mathrm{sin}\mathrm{\Omega }t+q_M(0)+\sqrt{3}m(1\mathrm{cos}\mathrm{\Omega }t)q_Y(0),`$
$`q_Y(t)`$ $`=`$ $`\mathrm{}\mathrm{\Delta }_{Y0}^2{\displaystyle \frac{1+3m^2}{2\sqrt{2}}}p_N(0)\mathrm{sin}\mathrm{\Omega }t\mathrm{}\mathrm{\Delta }_{Y0}^2\sqrt{3}mp_M(0)\mathrm{sin}\mathrm{\Omega }t+\mathrm{}\mathrm{\Delta }_{Y0}^2p_Y(0)\mathrm{sin}\mathrm{\Omega }t+q_Y(0)\mathrm{cos}\mathrm{\Omega }t,`$
$`p_+^{}(t)`$ $`=`$ $`{\displaystyle \frac{3(1m^2)+(1+3m^2)\mathrm{cos}\mathrm{\Omega }t}{4\sqrt{3}}}p_N(0)+{\displaystyle \frac{1+mm\mathrm{cos}\mathrm{\Omega }t}{\sqrt{2}}}p_M(0)+{\displaystyle \frac{\mathrm{cos}\mathrm{\Omega }t}{\sqrt{6}}}p_Y(0){\displaystyle \frac{\mathrm{sin}\mathrm{\Omega }t}{\sqrt{6}\mathrm{}\mathrm{\Delta }_{Y0}^2}}q_Y(0),`$
$`p_0^{}(t)`$ $`=`$ $`{\displaystyle \frac{3(1+m^2)(1+3m^2)\mathrm{cos}\mathrm{\Omega }t}{2\sqrt{3}}}p_N(0)\sqrt{2}m(1\mathrm{cos}\mathrm{\Omega }t)p_M(0){\displaystyle \frac{2\mathrm{cos}\mathrm{\Omega }t}{\sqrt{6}}}p_Y(0)+{\displaystyle \frac{2\mathrm{sin}\mathrm{\Omega }t}{\sqrt{6}\mathrm{}\mathrm{\Delta }_{Y0}^2}}q_Y(0),`$
$`p_{}^{}(t)`$ $`=`$ $`{\displaystyle \frac{3(1m^2)+(1+3m^2)\mathrm{cos}\mathrm{\Omega }t}{4\sqrt{3}}}p_N(0){\displaystyle \frac{1m+m\mathrm{cos}\mathrm{\Omega }t}{\sqrt{2}}}p_M(0)+{\displaystyle \frac{\mathrm{cos}\mathrm{\Omega }t}{\sqrt{6}}}p_Y(0){\displaystyle \frac{\mathrm{sin}\mathrm{\Omega }t}{\sqrt{6}\mathrm{}\mathrm{\Delta }_{Y0}^2}}q_Y(0),`$
$`q_+^{}(t)`$ $`=`$ $`\left[\sqrt{3}Nt\left({\displaystyle \frac{(1m)\epsilon }{1+m}}+4(c_0+c_2)𝒪_{\varphi \theta }\right)+{\displaystyle \frac{\sqrt{3}\mathrm{}\mathrm{\Delta }_{Y0}^2}{8}}(1+3m^2)(1+m)^2\mathrm{sin}\mathrm{\Omega }t\right]p_N(0)`$
$`+\left[{\displaystyle \frac{2\sqrt{2}\epsilon }{1+m}}Nt{\displaystyle \frac{3m(1m)^2}{2\sqrt{2}}}\mathrm{}\mathrm{\Delta }_{Y0}^2\mathrm{sin}\mathrm{\Omega }t\right]p_M(0)+\mathrm{}\mathrm{\Delta }_{Y0}^2{\displaystyle \frac{3(1m)^2\mathrm{sin}\mathrm{\Omega }t}{2\sqrt{6}}}p_Y(0)`$
$`+{\displaystyle \frac{1}{\sqrt{3}}}q_N(0)+{\displaystyle \frac{1}{\sqrt{2}}}q_M(0)+\left[{\displaystyle \frac{16m+3m^2}{2\sqrt{6}}}+{\displaystyle \frac{(1m)^2}{2\sqrt{2}}}\mathrm{cos}\mathrm{\Omega }t\right]q_Y(0),`$
$`q_0^{}(t)`$ $`=`$ $`\left[\sqrt{3}Nt\left({\displaystyle \frac{1+m^2}{1m^2}}\epsilon +4(c_0+c_2)𝒪_{\varphi \theta }\right)+{\displaystyle \frac{\sqrt{3}}{8}}\mathrm{}\mathrm{\Delta }_{Y0}^2(3m^42m^21)\mathrm{sin}\mathrm{\Omega }t\right]p_N(0)`$
$`+\left[{\displaystyle \frac{2\sqrt{2}m\epsilon }{1m^2}}Nt+\mathrm{}\mathrm{\Delta }_{Y0}^2{\displaystyle \frac{3m(1m^2)}{2\sqrt{2}}}\mathrm{sin}\mathrm{\Omega }t\right]p_M(0)+\mathrm{}\mathrm{\Delta }_{Y0}^2{\displaystyle \frac{3m^23}{2\sqrt{6}}}\mathrm{sin}\mathrm{\Omega }p_Y(0)`$
$`+{\displaystyle \frac{1}{\sqrt{3}}}q_N(0)+\left[{\displaystyle \frac{1+3m^2}{2\sqrt{6}}}+{\displaystyle \frac{3m^23}{2\sqrt{6}}}\mathrm{cos}\mathrm{\Omega }t\right]q_Y(0),`$
$`q_{}^{}(t)`$ $`=`$ $`\left[\sqrt{3}Nt\left({\displaystyle \frac{(1+m)\epsilon }{1m}}+4(c_0+c_2)𝒪_{\varphi \theta }\right)+{\displaystyle \frac{\sqrt{3}\mathrm{}\mathrm{\Delta }_{Y0}^2}{8}}(1+3m^2)(1+m^2)\mathrm{sin}\mathrm{\Omega }t\right]p_N(0)`$ (112)
$`\left[{\displaystyle \frac{2\sqrt{2}\epsilon }{1m}}Nt+{\displaystyle \frac{3m(1+m)^2}{2\sqrt{2}}}\mathrm{}\mathrm{\Delta }_{Y0}^2\mathrm{sin}\mathrm{\Omega }t\right]p_M(0)+\mathrm{}\mathrm{\Delta }_{Y0}^2{\displaystyle \frac{3(1+m)^2\mathrm{sin}\mathrm{\Omega }t}{2\sqrt{6}}}p_Y(0)`$
$`+{\displaystyle \frac{1}{\sqrt{3}}}q_N(0){\displaystyle \frac{1}{\sqrt{2}}}q_M(0)+\left[{\displaystyle \frac{1+6m+3m^2}{2\sqrt{6}}}+{\displaystyle \frac{(1+m)^2}{2\sqrt{2}}}\mathrm{cos}\mathrm{\Omega }t\right]q_Y(0).`$
## Appendix F The variance of macroscopic condensate spin in Cartesian coordinates

In Cartesian coordinates, the variances of the condensate spin direction is given by

$`\left[\delta f_x(t)\right]^2`$ $`=`$ $`{\displaystyle \frac{1}{1m^2}}\left[3p_N^2(0)2\sqrt{6}mp_N(0)p_M(0)+2m^2p_M^2(0)\right],`$
$`\left[\delta f_y(t)\right]^2`$ $`=`$ $`{\displaystyle \frac{1m^2}{2\mathrm{}^2}}[q_M^2(0)+2\sqrt{3}mq_M(0)q_Y(0)+3m^2q_Y^2(0)]`$
$`+{\displaystyle \frac{4N^2\epsilon ^2t^2}{\mathrm{}^2(1m^2)}}[3m^2p_N^2(0)2\sqrt{6}mp_N(0)p_M(0)+2p_M^2(0)]`$
$`+{\displaystyle \frac{\sqrt{2}N\epsilon t}{\mathrm{}^2}}\times [\sqrt{2}(q_M(0)p_M(0)+p_M(0)q_M(0))\sqrt{3}m(q_M(0)p_N(0)+p_N(0)q_M(0))`$
$`+\sqrt{6}m(q_Y(0)p_N(0)+p_N(0)q_Y(0))3m^2(q_Y(0)p_M(0)+p_M(0)q_Y(0))],`$
$`\left[\delta f_z(t)\right]^2`$ $`=`$ $`2p_M^2(0).`$ (113)
## Appendix G Coefficients in the zero mode Hamiltonian for anti-ferromagnetic interactions

The coefficients of zero mode Hamiltonian (97) for anti-ferromagnetic interactions (with $`=0`$)

$`a`$ $`=`$ $`{\displaystyle \frac{(1+3n_0)\epsilon }{6n_0(1n_0)}}+\left[6c_0{\displaystyle \frac{8c_2n_0}{3(1n_0)}}\right]𝒪_{\varphi \theta },`$
$`b`$ $`=`$ $`{\displaystyle \frac{\epsilon }{1n_0}}+4c_2𝒪_{\varphi \theta },`$
$`c`$ $`=`$ $`{\displaystyle \frac{\epsilon }{3n_0(1n_0)}}+{\displaystyle \frac{4c_2n_0}{3(1n_0)}}𝒪_{\varphi \theta },`$
$`\alpha `$ $`=`$ $`{\displaystyle \frac{2(3n_01)\epsilon }{3\sqrt{2}n_0(1n_0)}}+{\displaystyle \frac{16c_2n_0}{3\sqrt{2}(1n_0)}}𝒪_{\varphi \theta },`$
$`\eta `$ $`=`$ $`3c_2n_0(1n_0)𝒪_{\varphi \varphi }/\mathrm{}^2.`$ (114)


text:
# Potential Energy Landscape in Lennard-Jones binary mixture model

## I Introduction

In recent years a considerable research effort has been devoted to understand the complex phenomenology of supercooled glass forming liquids and in particular the enormous increase of relaxation times and viscosity by many order of magnitude upon decreasing the temperature. It was GoldsteinGoldstein69  who first related the behaviour of glass formers to the underlying Potential Energy Landscape (PEL) and proposed to characterize the dynamics of the system through the motion of a point in the complex high-dimensional PEL ($`3N`$-dimensional if $`N`$ is the total number of particles). Further, he suggested to focus onto the PEL local minima where the system is supposed to be trapped at low enough temperatures. At low temperatures, but above the glass transition, the phase space spanned by the system at equilibrium can ideally be characterized by two different types of processes: a ”fast” relaxation into local minima (or basins) and the ”slow” relaxation due to the crossing between loosely connected basins (the so called ”hopping” from basin to basin even if, in the classical mechanics PEL language, there exist only tortuous and crooked ways connecting two admissible basins). Goldstein focused the attention on the minima of PEL, but before the system is spending much of the time in a given basin, there are situations where ”fast” and ”slow” relaxations are not recognizable, diffusion is still high and a lot of minima are visited in a continuos way. In this case, it was searched to relate the diffusion dynamics to unstable directions in equilibrium (or instantaneous) PEL configurations, i.e. in the context of Instantaneous Normal Mode approachKeyes97 ; Donati00 , however other special points than PEL minima, are found to characterize better the dynamics of the system: they are the minima of the modulus square of the gradient of the potential energy, i.e. saddles (absolute minima) and what we call ”quasi-saddles” (local minima). Indeed the unstable directions of a saddle give some piece of information about the number of different minima are existing in that PEL region and therefore about the diffusion processes. In the past, the relevance of saddles for the dynamics of glassy systems has been already recognized in the context of mean-field spin-glassesKurchan96 ; Cavagna98 . A large number of recent works has renew interest in PEL and its characterization through saddle and minima points to better understand the continuos transition from liquid to glassPRL85-5356 ; PRL85-5360 ; Shah01 ; Cavagna01 ; Grigera02 .

The purpose of the present work is to investigate in detail the statistical properties of the PEL of a binary Lennard-Jones system from the liquid down to the supercooled regime near the structural arrest. To do that, a set of temperature above the glass transition have been considered and for each of them a set of independent configurations has been examined. Nearest saddles and minima have been searched and the normal mode analysis has been performed on all configurations.

## II Numerical computations

The system under consideration is a binary mixture of classical particles. The system is equal to that considered by Kob and AndersenKob95 , apart from a detail in the truncation of the potential. The mixture is made of $`256`$ particles, $`205`$ of $`A`$ and $`51`$ of $`B`$ (about $`80\%`$ of $`A`$ and $`20\%`$ of $`B`$). Both particles have the same mass and interact each other via a Lennard-Jones potential, i.e. $`V_{\alpha \beta }\left(r\right)=4ϵ_{\alpha \beta }[(\sigma _{\alpha \beta }\mathrm{}r)^{12}`$ $`(\sigma _{\alpha \beta }\mathrm{}r)^6]`$ with $`\alpha `$,$`\beta \{A,B\}`$. The parameters of the potentials are: $`ϵ_{AA}=1.0`$, $`\sigma _{AA}=1.0`$, $`ϵ_{BB}=0.5`$, $`\sigma _{BB}=0.88`$, $`ϵ_{AB}=1.5`$, $`\sigma _{AB}=0.8`$. In the following all the results will be given in reduced units, i.e. length in units of $`\sigma _{AA}`$, energy in units of $`ϵ_{AA}`$, temperature in $`k_B\mathrm{}ϵ_{AA}`$, and time in units of $`(m_A\sigma _{AA}^2\mathrm{}ϵ_{AA})^{1/2}`$.

The interaction potential at long distance is tapered between $`0.95r_2=r_1`$ $`rr_2=2.56\sigma _{AA}`$ with the following fifth-order smoothing function $`𝒯(r)=1+(r_1r)^3(6r^2+(3r+r_1)(r_15r_2)+10r_2^2)\mathrm{}(r_2r_1)^5`$. In this way the potential, the forces and their derivatives are continuos and it is possible to keep the energy constant better than $`1/10^5`$ over $`100`$ millions of molecular dynamics (MD) time steps near the critical temperature. The MD was performed in NVE ensemble using the leapfrog algorithm with a time step of $`\mathrm{1.5\; 10}^3`$ at temperature above $`T=1.0`$ and $`\mathrm{2\; 10}^3`$ at lower temperatures. A neighbor-list was used to speed up the calculation.

In their paper, Kob and AndersenKob95  truncated and shifted the potential at a cutoff distance of $`2.5\sigma _{\alpha \beta }`$. In the uniform density approximation, our tapering correspond to an average cutoff distance of about $`2.5\sigma _{AA}`$, and the resulting potential energy is shifted (downward) only by about 0.56 energy units with respect to the values of the quoted referencephishift .

First the configurations were produced at high temperature and then thermalized at the highest investigated temperature ($`T=2.0`$). The initial configurations at lower temperatures were obtained in sequence by cooling and thermalizing a set of different temperature ($`1.8`$, $`1.6`$, $`1.4`$, $`1.2`$, $`1.0`$, $`0.9`$, $`0.8`$, $`0.7`$, $`0.6`$, $`0.5`$, $`0.48`$, $`0.46`$). Not all the configurations were analyzed in full details.

For each temperature, $`1000`$ independent configurations were stored and then analyzed. The minimum number of steps between two subsequent configurations to be independent was determined in test runs. Two subsequent configurations are classified as independent if different potential minima are always obtained practically. We used a separation of $`1000`$ steps at high temperatures down to $`T=1.0`$, then we increased this number following the viscosity trend, so we used $`10^4`$ steps at $`T=0.7`$ and $`10^5`$ steps at $`T=0.46`$. The procedure is similar to that used to study the inherent dynamics by Schroeder et al.Schroeder00

To characterize the states of the supercooled binary mixture, we searched the nearest saddle point (”*W*”) and the nearest minimum (”*V*”) of the PEL from the initial equilibrium (”*I*”) configurations. The saddle (or quasi-saddle as explained after) configurations correspond to minimizing the sum of the squared forces, i.e. $`W=\left|V\right|^2`$ where $`V`$ is the sum of all interaction potential of our system and $``$ is the gradient on all the particle coordinates. If $`\left|V\right|^2`$ is really vanishing, i.e. is an absolute minimum, the configuration correspond to a real saddle, otherwise we call it quasi-saddle since usually there is only one direction with zero or nearly zero eigenvalue and non zero forcecm0108310 . The minima of $`V`$ are always absolute minima of $`\left|V\right|^2`$.

To obtain good nearest configurations of $`W`$ and $`V`$ minimum from an initial equilibrium configuration of our system of 256 particles is a very stiff problem. The details of the two search procedure we used, is too long and will be described elsewhereMS03 . Here we want to stress that various standard algorithms have been employed during the search and various criteria have been adopted for switching among the different algorithms. From the numerical point of view, the obtained minima are quite satisfactory. Indeed at the start, the value of $`\left|V\right|^2`$ is of the order of $`10^7`$ (in internal MD units) while at the minimum of $`W`$ is in the range $`10^1÷10^2`$ for quasi-saddles and in the range $`10^8÷10^2`$ for true-saddles. At the minimum of $`V`$ is always less than $`10^3`$ (in the range $`10^8÷10^2`$). In the worst case, the value of $`\left|V\right|^2`$ is reduced by a factor of $`10^5`$. This happens when the flex point is on a high sliding surface.

In figure 1 we plot the distribution of $`\left|V\right|^2`$ at the $`W`$ minima for two extreme values of temperature. It should be noted that: i) the temperature strongly affects the distribution; ii) the number of true-saddles (for which $`\left|V\right|^2`$ is vanishing) is strongly decreasing with temperature; iii) the force residuals, i.e. $`\left|V\right|^2`$ values, are larger at low temperatures in contrast with what one can expect intuitively. Moreover the temperature variation of the distribution is nearly confined to the low temperature side (below $`T=1.0`$). The second point is the reason why it becomes more and more difficult to analyze the low temperature data by using only these configurations which have the $`W`$ minima as true-saddles. Therefore we choose to analyze all the configurations and to inspect if true-saddles and quasi-saddles share the same pieces of information.

## III Analysis of the results

### III.1 Behaviour of potential energy at minima and saddles

In figure 2 we have reported the behaviour of the potential energy versus temperature of the quasi-saddles (”*W*”) and minima (”*V*”). As you can see, the potential energy of ”*W*” is decreasing rapidly with temperature while that of ”*V*” is practically constant, as already reported in ref. PRL85-5356 . Close and closer to the temperature of structural arrest, even the ”*V*” potential energy becomes a decreasing function of the temperature, as it has been already shown by Sastry et al.Sastry99

In the right panel of the same figure we have shown the Cartesian distances (in box units) of ”*V*” and ”*W*” from the starting equilibrium configurations. In that case the temperature behaviour of the two kind of minima is reversed between ”*V*” and ”*W*”, i.e. the average distances of quasi-saddles from equilibrium configurations are practically temperature independent while those of minima strongly decrease with decreasing temperature. This kind of behaviour can be expected because the potential energy of the starting equilibrium configurations is known to be nearly proportional to temperaturePRL85-5356  and configurations that are not very ”distant” from equilibrium are reasonably dependent on temperature.

In figure 3 we plot the number of negative eigenvalues of the Hessian matrix at quasi-saddle points versus temperature. This number goes fast to zero at the critical temperature, i.e. the temperature of structural arrest, as already found by a similar analysis performed in ref. PRL85-5356 ; PRL85-5360 . A solid vertical line in the figure represent the value of the mode coupling critical temperatureMCT92  (T$`{}_{c}{}^{}0.435`$) derived from the power law behaviour of the diffusion constantKob95 . The same value is obtained by a power law fitting of saddle order versus temperature (the solid line in figure 3) . This confirms the number of negative eigenvalues of quasi-saddles (or saddles) to be strictly connected with the diffusion constant. In contrast, a negative eigenvalue of an equilibrium configuration, i.e. a negative curvature of the instantaneous potential energy, is not always connected with a path joining two different minima, i.e. a diffusion path.

### III.2 Quasi-saddles characteristics

As already noted in ref. PRL85-5360 , the quasi-saddles at different temperatures explore different regions of the PEL. Therefore at low temperatures we have low values of the potential energy of ”*W*” together with low numbers of negative eigenvalues, as shown in figure 4.

We underline that the potential energy itself is not sufficient to determine the number of negative eigenvalues of the Hessian matrix. Indeed even at the same value of the potential energy the number of negative eigenvalues is substantially different for quasi-saddles (”*W*”) and initial equilibrium configurations (”*I*”), as shown in figure 5 where we plotted the potential energy of ”*I*” and ”*W*” versus the number of negative eigenvalues at different temperatures. From the values and distribution in energy of ”*I*” and ”*W*” at the same temperature, we can easily derive that the enormous decrease of diffusion (and increase in viscosity) lowering the temperature in supercooled regime, is due essentially to entropic barriers, i.e. to difficulties in reaching the saddle points, and not to activated processes Sampoli03 .

### III.3 Spectral densities of potential minima, saddles and initial configurations

To better understand the modifications of the PEL around ”*I*”,”*W*”,”*V*” configurations that occur lowering the temperature we have study the frequency distribution of the vibrational states derived from the Hessian matrix of the potential energy. To visualize all the spectral features in a single plot, the eigenfrequencies of unstable directions (i.e. the square roots of negative eigenvalues) are reported as negative frequencies.

In figure 6 we reported the histogram of the average vibrational density of states (DOS) at two extreme temperatures. From the figure we can infer that by lowering the temperature all the three DOS will collapse into a single one (when the glass will be confined to a single minimum). In the supercooled regime, approaching Tc, the ”*W*” DOS from quasi-saddles becomes more and more similar to the ”*V*” DOS from potential minima. At T$`=0.46`$, the ”*W*” DOS is nearly coincident with the ”*V*” DOS while the ”*I*” DOS of initial equilibrium configuration still remain separated. We want to underline that in the channel of the histogram around zero frequency, we found always the three eigenvalues connected to the three translational invariants; for quasi-saddles (true-saddles apart) another zero frequency eigenvalue is found due to the change of curvature in one direction of the potential energy hypersurface. This is the reason why the zero channel of the ”*W*” DOS has nearly four eigenvalues.

In figure 7, the three DOS calculated from only the true-saddle configurations, i.e. those configurations for which the forces at ”*W*” points are really vanishing, are shown and compared with the DOS averaged over all configurations. We want to stress that within the statistical uncertainties the three DOS are not dependent on the chosen subset. This is not obvious for ”W” DOS and it is an important result of our study: at least with respect to the distribution of vibrational frequency, true-saddles and quasi-saddles share the same information.

Even if the ”*W*” and ”*V*” DOS seem to be very similar near Tc (see lower panel of figure 6) there are yet significant differences at low frequency. In figure 8, we reported the reduced density of states (RDOS) at low frequencies for the three configuration at T$`=0.48,`$ i.e. the DOS divide by the frequency square. In the Debey solid approximation the reduced density is a constant, while the well-known ”Boson peak” is found in real glasses. In the supercooled regime of figure 8, a peak is clearly visible in the ”*V*” RDOS, starts to appear in ”*W*” RDOS, while is (yet) completely absent in the ”*I*” RDOS.

## IV Conclusions

We have investigated the dynamics of a model glass-forming liquids in terms of its potential energy landscape by investigating numerous independent configurations at different temperatures in the liquid and supercooled regime. Potential energy minima, quasi-saddles and true-saddles are derived from equilibrium configurations. The temperature behavior of potential energy, saddle order, Cartesian distance, and vibrational density of states for equilibrium, quasi-saddle and minimum configurations has been studied. In the average the quasi-saddles lie on a more sliding hypersurface at low than high temperatures. The number of true saddles is a decreasing function of temperature, so it is important to establish that it is possible to extract also information about the dynamics of the system from the quasi-saddles configurations. It has been demonstrated that the vibrational density of states is practically the same for true- and quasi-saddles. At different temperatures, different regions of the PEL are statistically visited by the system, due to the fact the accessible configuration volume is exponentially growing with the temperature. The potential energy and the number of unstable directions (i.e. negative eigenvalues) are statistically correlated both for saddles and equilibrium configurations, but at the same potential energy saddles and equilibrium configuration have different number of negative eigenvalues, at least in the investigated liquid-supercooled region. An excess density of states at low frequency is always present in the investigated system, but a ”Boson peak” like structure is clearly visible only in the DOS derived from potential energy minima.

This work was realized with the financial support of MURST.