Fragmentation, domain formation and atom number fluctuations of a two-species Bose-Einstein condensate in an optical lattice

We theoretically study the loading of a two-species Bose-Einstein condensate to an optical lattice in a tightly-confined one-dimensional trap. Due to quantum fluctuations the relative inter and intra species phase coherence between the atoms and the on-site atom number fluctuations are reduced in the miscible regime. For the immiscible case the fluctuations are enhanced and the atoms form metastable interleaved spatially separated domains where the domain length and its fluctuations are affected by quantum fluctuations.

In this paper we study both numerically and analytically the effects of quantum and thermal fluctuations on a two-species BEC when the condensates are confined in an optical lattice in a highly-elongated 1D trap and the lattice potential is slowly turned up. In singlespecies bosonic atomic gases the interplay between enhanced quantum fluctuations in an optical lattice and the repulsive inter-atomic interactions has experimentally been shown to result in strongly reduced atom number fluctuations and the loss of phase coherence between the atoms in different lattice sites [25,26,27,28,29,30,31,32,33]. The reduction in atom number fluctuations has been exploited in preparation of spin-squeezed states [29] that are suitable for quantum-enhanced interferometry. The phase separation dynamics of a twospecies harmonically-trapped BEC was experimentally observed in two immiscible hyperfine components of 87 Rb [1], due to the long lifetimes of the two-fluid system that results from a fortuitous cancellation of the scattering lengths [34]. Controllable spatial separation dynamics has recently been observed in a mixture of 85 Rb and 87 Rb atomic BECs by tuning the interspecies interactions with a magnetic Feshbach resonance [2] and by introducing a linear electromagnetic coupling between the two internal states [35] that creates an effective dressed state description for the two components [36,37].
In our simulations the atoms are initially confined in a shallow lattice and we continuously turn up the lattice potential. The enhanced effective interactions result from the reduced hopping amplitude of the atoms along the lattice and quantum fluctuations become more dominant in a deep lattice. We calculate the on-site atom number fluctuations in individual lattice sites and the inter and intra species relative phase coherence between the atoms in different sites. The numerical results are compared with the analytically calculated values that we derive in the appendix using the Bogoliubov theory. Even in the miscible regime of the two-species BEC system quantum fluctuations eventually destroy the longrange coherence of the atoms along the lattice, fragmenting the condensates. The inter and intra-species relative phase coherence between the atoms even in the adjacent sites is notably reduced in deep lattices, but the inter-species coherence remains higher close to the onset of the phase separation instability. We find that the repulsive inter-species interactions increase the inter-species relative phase coherence, but have only a weak effect on the intra-species coherence and the on-site atom number fluctuations. The coherence typically stabilizes to a non-vanishing finite value after the lattice ramping and we evaluate its spatial dependence along the lattice, demonstrating a clearly reduced spatial coherence length of the system. In the dynamically unstable regime we find considerably enhanced atom number fluctuations, stronger loss of phase coherence, and the spontaneous formation of metastable configurations of interleaved domains of the two spatially-separated components. We calculate quantum mechanical expectation values and uncertainties of the domain length and find that they depend on the strength of quantum fluctuations, deviating from the classical mean-field values.
The experimentally observed phase separation dynamics of Ref. [2] in the uniform space was theoretically studied using the classical mean-field theory in Ref. [38]. A multiorbital wavefunction analysis was employed in Ref. [39] to demonstrate that stronger inter-species interactions lead to a shorter domain length. There has been an increasing interest in the experimental studies of spontaneous symmetry breaking and pattern formation also in other ultra-cold atomic systems [40,41] and the two-species condensate with a coupling between the two internal states [35] has been proposed as a system to study the Kibble-Zurek defect formation mechanism in phase transitions [42,43].
In order to perform efficiently the numerical simulations in a lattice, we develop an approximate method to describe the non-equilibrium dynamics of a two-species condensate mixture that is based on the truncated Wigner approximation (TWA) [44,45,46,47,48,49,50] of the stochastic phase space dynamics. The two-species BEC equilibrium state is solved within the classical Bogoliubov approximation where the excitations are evaluated in the tight-binding approximation of the uniform two-species lattice Hamiltonian. The amplitudes of the Bogoliubov phonon modes are then stochastically sampled according to a probability distribution given by the Wigner distributions of the ideal harmonic oscillators, as in the single-component TWA approaches [47]. Each stochastic realization of the initial state is propagated in time according to the classical mean-field dynamics, so that individual stochastic trajectories represent potential outcomes of single experimental runs and quantum mechanical expectation values and fluctuations are calculated from the ensemble averages of the stochastic dynamics. One of the advantages of the approximate two-species model is its simplicity and the possibility to calculate analytic solutions to the initial state mode functions. The analytic approach to linearized excitations can also be used to calculate approximate ground state properties of the two-species system. We use this in Appendix to evaluate the intra-species relative phase coherence and the on-site atom number fluctuations.

Classical mean-field equation
We assume that a two-species BEC is in a tightly-confined highly-elongated 1D trap, so that any density fluctuations along the radial direction perpendicular to the trap axis can approximately be ignored. Along the axial direction the atoms experience an optical lattice potential that is deep enough so that the atoms can be described in the tight-binding approximation in which case one trap mode per lattice site is included in the dynamics. The classical mean-field model then follows our previous descriptions [8,16] and the equation  that governs the dynamics of two component BEC is the two-component discrete nonlinear  Schrödinger equation (TCDNLSE), where J j (J j > 0) and ψ ( j) n denote the nearest-neighbour hopping amplitudes and the wavefunction amplitude at the lattice site n of the atoms of species j ( j = 1, 2), respectively. The nonlinearities are given by the interaction coefficients χ jk ( j, k = 1, 2) which are proportional to the onsite atom-atom interaction strengths and to the overlap of the lowest vibrational state wave functions φ ( j) n (the Wannier functions) of the two species in a given lattice site, i.e., and for j k. The inter and intra species scattering lengths are denoted by a jk and a j j , respectively, and N j is the atom number of the species j. Here the reduced atomic mass µ is given in terms of the atomic mass of the j th component m j as Since the number of atoms in each species is a conserved quantity, in the following we use the normalization where L denotes the number of lattice sites.

Parameter regimes of the simulations
We study the quantum dynamics of a two-species BEC in an optical lattice. The initial state of the TWA simulations is generated by calculating the classical Bogoliubov modes whose amplitudes are sampled stochastically. The time-evolution for each stochastic realization then follows TCDNLSE of Eq. (1). The simulations involve a large parameter space. There are three nonlinearities χ 11 , χ 22 , and χ 12 in Eq. (1) and the hopping amplitudes J k for the two species may differ. In addition, the two components can be moving with the different carrier momenta p k . The stochastic initial state fixes the atom numbers N 1 and N 2 (for given χ 11 , χ 22 ) and we can also vary the number of lattice sites L. The lattice potentials of the two components may be shifted with respect to each other and the two species could also experience different radial confinements. In the finite temperature examples we also vary the initial temperature T . In the following we will demonstrate how a simple analytic description for the initial state of the TWA simulations in terms of the Bogoliubov modes may be obtained whenever the two atom currents are equal, i.e., for J 1 sin(p 1 ) = J 2 sin(p 2 ). In order to demonstrate some basic effects of the two-species quantum dynamics in a lattice, we concentrate on a simple set of parameter values for which the mode functions have especially compact analytic expressions. In all the numerical simulations we consider condensates with zero centre-of-mass momenta p 1 = p 2 = 0 and We show that the particular choice for the set of parameters is by no means necessary, however, and that the general formalism with the analytic initial state derivation is more general. For the selected parameter set we may investigate the main physical phenomena of the two-species lattice dynamics: the effect of the phase separation dynamics can be controlled by the ratio χ 12 /χ 11 , quantum fluctuations and nonlinearity by χ 11 /N 2 1 and χ 11 /J 1 .

Classical Bogoliubov theory and stability analysis
We can find steady-state solutions for the TCDNLSE (1) that represent propagating plane waves [8,16] where ω j is given by The carrier wave momenta, p j are quantized according to p j = P j 2π L where P j is an integer that takes value in the interval [− L 2 , L 2 ). The linear stability analysis of the steady-state solution was performed in Ref. [8] using the classical Bogoliubov expansion. In the Bogoliubov approach the wavefunctions for each component in Eq. (1) are written as In the limit of weak perturbations, the system of equations for u ( j) q and v ( j) q may be expressed as an eigenvalue problem where σ z denotes the 2 × 2 Pauli spin matrix. The elements of the 4 × 4 matrix M q are obtained from the Bogoliubov linearization procedure [8]. The quasimomenta q = 2πQ L may be defined such that Q takes integer values except zero in the range [− L 2 , L 2 ). The eigenvalues of M q correspond to the normal mode (excitation) frequencies Ω q of the system that have simple analytic expressions when two BECs have the same atomic currents [8](J 1 sin(p 1 ) = J 2 sin(p 2 )). In that case we obtain [8,16] where represents a Doppler shift term of the excitation frequencies due to the superfluid flow, are defined in terms of the single-condensate normal mode frequencies ν j,q and is the spectrum of an ideal, non-moving BEC. The flow is stable if the frequencies in Eq. (11) are real for all q 0, otherwise there are small excitations in the system that grow exponentially in time. In these equations the normal mode frequencies of the two condensate species are coupled by the inter-species interactions and in the absence of the inter-species term, χ 12 = 0, we have two decoupled condensate spectra Eq. (20). The simplest case is obtained when both BECs are in the normal dispersion regime with p 1 , p 2 < π/2 and χ jk > 0.
In that case the instability condition for the modes q reads [8] where j k = J k cos p k . In the normal dispersion regime the instability first sets in for the modes for which |q| is small and therefore q 2 L ∝ 1/L. In the limit of a large lattice L → ∞, we then obtain the criterion for instability This stability condition is notably altered if one of the BECs exhibits anomalous dispersion due to superfluid flow [8]. If the interaction strengths in Eq. (17) are tuned in such a way that the instability is characterized by a single unstable mode, the two-component system can also be found in a state that is no longer dynamically stable but does not undergo a phase separation [16]. Instead, the two-species mixture exhibits a periodically appearing and disappearing vector soliton structure.
In this work we only consider initially stationary BECs with the vanishing condensate momenta p 1 = p 2 = 0. We also assume that the hopping amplitudes for the two BECs are equal J 1 = J 2 = J and ǫ 1,q = ǫ 2,q ≡ ǫ q = 4J sin 2 (q/2). This simplifies the stability analysis.
The eigenvalue system for the linear stability analysis in Eq. (10) can then be expressed as M q given by with the definition η j,q = U j j + ǫ q . The system exhibits two physical normal mode frequencies For a positive definite M q these are real indicating dynamical stability of the system. The corresponding dynamically stable eigenvectors ξ q satisfy the normalization condition ξ † q σξ q = 1. The eigenvalues −Ω ± q of Eq. (10) represent unphysical solutions with the corresponding eigenvectors satisfying the negative normalization ξ † q σξ q = −1. The BEC system becomes dynamically unstable when the normal mode frequencies exhibit nonvanishing imaginary parts, indicating perturbations that grow exponentially in time. The corresponding eigenvectors satisfy ξ † q σξ q = 0. The rate at which the instability sets in depends on the magnitude of the imaginary part of the eigenfrequency.
We can solve the eigenvectors of Eq. (10) analytically. The expressions for the mode functions notably simplify when we consider the case χ 11 = χ 22 . We then obtain (for χ 12 0) and u (2) q,± = ±u (1) q,± , v (2) q,± = ±v (1) q,± . Here we have assumed, for the notational simplicity, that ∆ 12 > 0. In TWA simulation we generate the initial state noise for the configuration that is dynamically stable (we specifically consider thermal equilibrium states). For such states the normal mode frequencies are all real and the eigenmodes satisfy the normalization condition The two-species BEC normal modes describe the dynamics of mixing between the two components as well as excitations of the total density in the system. In nonlinear regime we have χ j j ≫ J. If we also have χ 2 12 ≃ χ 11 χ 22 , one of the frequencies approaches zero corresponding to the phase separation instability Eq. (18). We obtain in that case Ω 2 q,+ ≃ ν 2 1,q + ν 2 2,q and Ω 2 q,− ≪ ν 2 1,q , ν 2 2,q . The low energy excitations then correspond to the mixing of the two species with only a weak variation in the total density of the two-species condensate.
In both stable and unstable regimes of the two-species mixture we can investigate the degree of overlap between the two species. We define the overlap integral of the wavefunctions for the mixture as In the stable regime Eq. (25) describes the spin excitations of the two-component system. When the two-species interaction strengths satisfy the condition Eq. (17), the two-species system is dynamically unstable and undergoes phase separation, resulting in strongly reduced overlap integral values. A measure that can particularly well identify the phase-separation of the densities of the two species may be calculated from the sum

Truncated Wigner approximation
In the TWA simulations we calculate ensemble averages of stochastic trajectories for which the time evolution follows the classical mean-field theory, but in each realization the initial state is stochastically sampled from a Wigner distribution that approximately synthesizes the quantum statistical correlations of the initial state. Approaches introduced in the TWA initial state generation in single-component BECs involve evaluating the initial state correlations within the Bogoliubov approximation [51,47] or by solving the ground state and the excited state populations self-consistently within the Hartree-Fock-Bogoliubov approximation [33]. In order to implement the TWA phase-space model in a two-component BEC system we similarly assume that the two component stochastic fieldsψ ( j) obeys the classical field equations similar to the TCDNLSE (Eq. 1) For the stochastic initial state generation we introduce an approximate model based on the classical Bogoliubov theory, described in the previous section, that provides simple analytically solvable mode functions. We write the both components j = 1, 2 as N jψ where φ ( j) n denotes the normalized ground state solution of the BEC component j and the excited-state fluctuations are given bȳ Here the eigenmodes (q, ±) correspond to the eigenfrequencies Ω ± q of Eq. (20). The mode amplitudes α ( j) 0 , α ( j) q,η are stochastically sampled from the Wigner distribution of harmonic oscillators as explained below. We consider a two-species system with equal populations N 1 = N 2 ≡ N, with the interaction strengths satisfying χ 11 = χ 22 . In that case the modes u ( j) q,η and v ( j) q,η are given by Eqs. (22) and (23).
Using the field decomposition Eq. (28) and the mode functions, Eqs. (22) and (23), we then stochastically sample the amplitudes of the mode functions by treating them as ideal harmonic oscillators whose distributions are determined by the corresponding Gaussian Wigner function [52] where ξ ± q ≡ Ω ± q /2k B T . The stochastic mode function amplitudes α ( j) q,ν produce the ensemble averages wheren is the usual Bose-Einstein distribution function. The factor 1/2 in Eq. (31) results from the Wigner distribution that returns symmetrically ordered expectation values, providing the vacuum noise in each mode. For each stochastic realization the total number of excited-state atoms varies according to with the average number given by The ground state amplitudes α ( j) 0 fluctuate in each stochastic realization [50,53]. The groundstate atom number is then obtained from the fixed total atom number N in each atomic species,

Turning up the optical lattice
We solve the stochastic dynamics during the turning up of the lattice potential for given initial conditions. We study the response of the system to the ramping so that both the tunneling coefficients and nonlinearity are time-dependent. We assume that at all times J 1 = J 2 = J and χ 11 = χ 22 . For simplicity, we consider a situation where the both species have an equal mass m, so that the recoil frequencies ω R are equal and are given by where d denotes the lattice spacing. In a deep lattice the hopping amplitude J is then approximately given by [54] Here s denotes the lattice height in the units of the lattice photon recoil energy. In a tightlyconfined elongated 1D trap the atoms are assumed to be confined in the radial vibrational ground state. If the radial confinement is the same for the both species, we obtain where Ω ⊥ denotes the trapping frequency of the radial confinement and l s = ( /mΩ s ) 1/2 , where Ω s ≃ 2s 1/2 ω R is the axial trap frequency at the lattice site minimum. In our simulations we numerically solve the time-evolution using the split-step method [55]. As the lattice is turned up the hopping amplitude J rapidly decreases according to Eq. (36) and the interaction strength χ i j slowly increases according to Eq. (37), due to the stronger confinement of atoms in individual sites. We increase the lattice height linearly at the rate δ, so that the lattice height satisfies s(t) = s i + δt, where s i denotes the initial height. We choose δ = 2 × 10 −3 ω R and s i = 2, resulting in the initial value of J ≃ 0.22ω R . Unless it is stated explicitly, the length of the lattice is L = 64 and the number of atoms in each species in each run is taken to be N/L = 40. In most cases we choose the initial value for the intraspecies interaction strengths The initial state fluctuations are evaluated within the Bogoliubov approximation and the interaction strengths are selected in such a way that the excited state population remains low. For χ 11 = χ 22 ≃ 0.60ω R and χ 12 /χ 11 = 0.1 at T = 0, the number of atoms initially depleted from the ground state due to quantum fluctuations N (1) e = N (2) e ≃ 70 (corresponding to 2.7% depleted fraction). All the presented simulation results are for T = 0, except the finite temperature cases studied in Sec. 3.2.4.
We study the effect of ramping of the lattice both in the stable (χ 12 √ χ 11 χ 22 ) and unstable (χ 12 √ χ 11 χ 22 ) regime and write the inter-species interaction strength χ i j (for i j) as The nonlinearity corresponding to the inter-species interactions is tuned by varying the parameter γ in Eq. (38). Here γ 1 corresponds to the stable regime while γ 1 implies the dynamical phase-separation instability. In cases where we study the two-species system in the unstable phase separation regime, we use the initial two-species mixture that is dynamically stable, but change the value of γ from the stable to unstable regime immediately after the lattice ramp. For a fixed nonlinearity quantum fluctuations are enhanced by reducing the atom number (and correspondingly increasing the scattering length). For large atom numbers and in shallow optical lattices quantum effects are weak and the system can be accurately described by classical mean-field theory. As the lattice is turned up, the effective interactions become stronger and quantum fluctuations are enhanced.
The ramping of the lattice is adiabatic if the rate of change of the parameters in the Hamiltonian is slow compared to the lowest collective excitation frequency [56,47]. The fastest varying parameter during the ramping is the hopping amplitude and we require for an adiabatic ramp If this condition is not satisfied, the system can be excited from its ground state during the turning up process of the lattice.

Validity of the parameter regimes
Our chosen set of parameters, explained in Sec. 2.2, captures the essential features of the condensate fragmentation, reduced atom number fluctuations and the domain formation. For the phase separation dynamics, the important condition is that of the dynamical instability χ 2 12 χ 11 χ 22 and the precise ratio χ 11 /χ 22 is less relevant. Experimentally, the interaction strengths χ i j can be controlled using two-species Feshbach resonances [2] or by introducing a linear electromagnetic coupling between the two internal states [35]. The interaction parameter χ i j incorporates the atom numbers and the ratio χ 11 /χ 22 may also be tuned by changing the relative atom population of the two condensate components. The intra-species interaction strength χ 12 can be controlled in a spin-dependent optical lattice [20] by changing the relative lattice positions of the two species and therefore modifying the spatial overlap integral between the lattice site wavefunctions.
The two-species condensate system may be realized by using two different hyperfine levels of the same atom [1,2] or by trapping two entirely different atoms, e.g., a 41 K-87 Rb mixture [21]. In the case of far-detuned optical lattice, the potential experienced by the atoms in two different hyperfine levels of the same atom is typically the same, resulting in the identical values for the hopping amplitudes.
The atom dynamics can be described by a 1D model if the frequency of the radial trapping potential Ω ⊥ is larger than the chemical potential of the atoms ω j in Eq. (8) and the thermal energy k B T . The typical values used in the numerics are χ 11 = χ 22 ≃ 0.60ω R at s = 2, corresponding to χ 11 ≃ 1.3ω R at s = 40 and Ω ⊥ a 11 /d ≃ 7.9 × 10 −5 ω R . We therefore obtain for the requirement of the 1D dynamic description Ω ⊥ 1.3ω R /L and d/a 11 250. The first expression yields in terms of the lattice spacing d 0.3l ⊥ where l ⊥ = √ /(mΩ ⊥ ) denotes the radial width of the lattice site mode function. For the scattering length a 11 ≃ 5nm we obtain from the second inequality d 1.3µm. In 1D lattices the lattice spacing can be easily controlled by adjusting the angle between the counter-propagating lasers that form the standing-wave pattern, so that d = π/[k sin(θ/2)] for the intersection angle θ and wavenumber k. In the simulations we assume 40 atoms per site and we require that the atom density is sufficiently low so that the inelastic atom losses remain weak. The three-body loss rate of the atoms at the site l may be approximated by where K 3 denotes the three-body recombination rate. In order to have a weak three-body loss rate n 3 l Γ/ω R ≪ 1, the radial trap frequency should be sufficiently weak. The parameters depend on the particular atom and the hyperfine state. For instance, for 87 Rb in the |F = 1, m F = −1 state we have K 3 ≃ 5.8 × 10 −30 cm 6 /s [57] and we obtain for the condition of the three-body loss rate to be weak Ω ⊥ ≪ 1.2 × 10 5 /s. This condition can be satisfied when the system dynamics is strictly 1D. It should be also noted, however, that even in elongated traps that are not tightly confined, 1D numerical model can provide a good qualitative description of the atom dynamics in the lattice [33].
In the tight-binding approximation we assume that only the lowest energy band is occupied and one mode function per lattice site is sufficient to represent the dynamics. This approximation becomes better in deeper lattices and provides a reasonable description for s 2 [54]. We also require that the nonlinearity is smaller than the energy gap between the lowest two energy bands ∼ Ω s = 2 √ s ω R , so that the higher bands are not occupied. This yields 2 √ sω R χ i j /L which is well satisfied for all the studied lattice heights. Experimentally, the atoms are trapped in a combined optical lattice and a harmonic trap. The harmonic potential introduces a non-uniform atom density. This influences the phase separation dynamics; the condensate component with a weaker nonlinearity energetically favours higher density regions close to the centre of the trap [1]. The phase coherence and number fluctuations in a harmonic trap depend on the spatial location with quantum and thermal fluctuations stronger close to the edge of the atom cloud [33]. Some other possible effects on the phase coherence are addressed in Sec. 3.2.1. A lattice with a uniform density and periodic boundary conditions may be realized in a toroidal trap with an optical lattice formed by the interference of two counter-rotating Laguerre-Gaussian laser beams [58].

Dynamically stable regime
We first consider the two-species BEC dynamics in an optical lattice in the dynamically stable regime of the spatially overlapping condensate mixture for γ < 1. This corresponds to the situation where the inter-species interaction is not strong enough to cause the phase separation of the two components and all the normal mode frequencies in Eq. (11) are real.

Condensate fragmentation and phase coherence
The atoms are initially confined in a shallow lattice and we continuously turn up the lattice potential. The effect of quantum fluctuations on atom dynamics in the lattice can be studied by calculating the phase coherence between the atoms in different lattice sites. When the lattice is turned up the hopping amplitude of atoms between adjacent sites rapidly decreases, resulting in the reduction of kinetic energy of the atoms and hence stronger effective interactions. Quantum fluctuations in the system are enhanced and the phase coherence between the atoms in different sites is destroyed as the condensates undergo fragmentation. We evaluate the loss of phase coherence by calculating the absolute value of the normalized relative intra-species phase coherence between the atoms in different sites, between the atoms of the same species j in sites k and l. Since we choose N 1 = N 2 = N and χ 11 = χ 22 , the two species have identical nonlinear properties that are spontaneously broken only due to nonlinear interactions, e.g., in the phase separation. We describe the relative phase coherence between the atoms in the two different species by In Fig. 1 we show the relative intra-species phase coherence between the atoms in two adjacent sites C 1 during and after the turning up of the lattice potential. The different curves correspond to different values of the final lattice height s = 10, 20, 30 and 40. The displayed cases have the inter-species interaction strength γ = 0.1 and γ = 0.5, defined by Eq. (38). We also show the corresponding stationary (averaged) values of the coherence C 1 that are obtained after the turning up of the lattice potential. These demonstrate how the intra-species coherence rapidly decreases as the lattice becomes deeper, indicating an increasing degree of fragmentation of the BEC as a function of the final lattice height. We only find a very weak dependence of the intra-species relative phase coherence C 1 on the inter-species interaction strength γ for the values of γ in the stable regime we considered in the simulations (from . The coherence C 10 decreases notably more rapidly than the coherence C 1 between the atoms in adjacent sites as a function of the lattice height. As in the case of C 1 , we find that C 10 depends only weakly on γ in the stable regime for the values we studied from γ = 0.1 to 0.75. γ = 0.1 to 0.95). In Ref. [21] the presence of 41 K was found to lead to lower visibility of the interference fringes of the 87 Rb in the two-species mixture in an optical lattice. In the experiment, however, the interactions increased the atom density of 87 Rb close to the trap centre due to the inhomogeneous trapping potential and drove the system closer to the onset of the Mott insulator transition, as also demonstrated in Ref. [59] where the adiabaticity and the visibility of interference fringes in loading of a bosonic mixture to an optical lattice was studied using Gutzwiller mean-field method. The intra-species long-range spatial coherence is shown in Fig. 2. We display the relative phase coherence C 10 between the atoms in one of the sites and in its 10th nearest neighbour site. This decays notably faster than the coherence C 1 between the atoms in the adjacent sites in Fig. 1. The spatial dependence of the relative phase coherence along the lattice is also shown in Fig. 2. We calculated the stationary, averaged values of C j for different j when the value of the coherence was stabilized after the end of the lattice ramp. The graphs show the decay of the spatial coherence along the lattice. The coherence exhibits a very slow decay for large values of j and remains high for the case of shallow lattices.
The numerically calculated values of the intra-species relative phase coherence may be compared to the analytic estimates obtained for the ground state of the optical lattice system in Appendix. The results for the nearest-neighbour phase coherence for the linearized Bogoliubov theory of the fluctuations in the ground state, displayed in Fig. A1 in Appendix, provide a good agreement with those obtained in the numerical TWA simulations, shown in Fig. 1. The long-range coherence along the lattice in the TWA numerics, however, decays more rapidly as a function of the spatial separation than in the case of the ground-state calculation.  The inter-species coherence C (12) 1 is shown in Fig. 3 for different final lattice heights and for different values of inter-species interactions, γ = 0.1, 0.25, 0.5 and 0.75. Although the intra-species coherence C 1 is not strongly affected by the inter-species interaction strength γ, the inter-species coherence C (12) 1 is very sensitive to γ even when the two-species mixture is miscible. In particular, the relative inter-species coherence C (12) 1 is enhanced due to the interspecies interactions when χ 12 is increased. C (12) 1 becomes high as the system approaches to the onset of the phase separation instability, as shown in stationary averaged values of Fig. 4 that are calculated after the turning up of the lattice. This is because the effective interactions in a perfectly overlapping two-species mixture are almost completely canceled out immediately below the onset of the instability for χ 11 ≃ χ 22 ≃ χ 12 . We also show the decay of the spatial coherence along the lattice by displaying C (12) j for different values of j in Fig. 4 (on right). The stationary, averaged values for the coherence are calculated after the end of the ramp when the coherence has stabilized.

Atom number fluctuations
for each species-j at the lattice site l. It is useful to scale the atom number fluctuations to those obtained in the Poissonian limit that correspond to the fluctuations resulting in an instantaneous splitting or in the splitting of a non-interacting gas so that the values n ( j) sqz < 1 indicate reduced on-site atom number fluctuations. The onsite and the relative atom number fluctuations between the atoms in different lattice sites were calculated within TWA for a single-species BEC and compared with experimental observations in Ref. [33] providing a good qualitative agreement.
In Fig. 5 we show the scaled on-site atom number fluctuations n sqz in one of the sites for the different interaction strengths. The stationary averaged values of atom number fluctuations are shown in Fig. 6. We show both the dependence of the fluctuations on the lattice height as well as on the inter-species interactions. The atom number fluctuations are strongly reduced as the final lattice height is increased. Strong suppression on number fluctuations for deep lattices is associated with the enhanced phase fluctuations and the condensate fragmentation, and the reduced atom number fluctuations correlate with the previously calculated values of the loss of phase coherence. We also find that the inter-species interactions generally enhance the atom number fluctuations.  The non-equilibrium dynamics of the TWA simulations may again be compared to the analytic estimates obtained for the ground state of the optical lattice system in Appendix. The agreement between the linearized ground-state results of Fig. A1 in Appendix and the TWA results of Figs. 5 and 6 is very good, indicating that the effects of nonadiabaticity in the lattice ramping on the on-site atom number fluctuations are weak. Fig. 7  the lattice potential results in a decrease in the lowest mode population. The effect is stronger when the lattice becomes deeper, indicating breakdown of the adiabaticity in turning up of the lattice. Consequently, the deeper the lattice, the larger the depletion in the lowest mode. The breakdown of the adiabaticity in deep lattices may be understood from Eq. (39) since the frequency of the lowest phonon mode is reduced as the lattice becomes deeper and it becomes progressively more difficult to turn up the lattice adiabatically. We may estimate the adiabaticity of the turning up of the lattice potential using the expression Eq. (39). In the TWA simulations the lattice is turned up at the rate s(t) = s i + δt, with s i = 2. For γ = 0.1 and δ = 2.0 × 10 −3 ω R we obtain Ω min q (t)/ζ(t) ≃ 10 at s = 2 and 0.7 at the end of the deep lattice ramp s = 40. The condition Ω min q (t)/ζ(t) ≃ 1 is reached at about s ≃ 34. The adiabaticity condition is much easier to violate close to the onset of the phase separation. For γ = 0.95 and δ = 2.0 × 10 −3 ω R we obtain Ω min q (t)/ζ(t) ≃ 4 at s = 2, reducing to 1 at s ≃ 12, and to 0.2 at s = 40.

Adiabaticity and excitations of relative atom populations In
In Fig. 8 we show the the relative phase coherence between the atoms in the adjacent lattice sites and the on-site atom number fluctuations for three different speeds of the lattice ramp, representing δ = 2.0 ×10 −3 ω R , 1.0 ×10 −3 ω R and 0.67 ×10 −3 ω R . For γ = 0.1 and the two slowest ramp cases the condition Ω min q (t)/ζ(t) ≃ 1 is never reached during the turning up of the lattice potential. For 0.67 × 10 −3 ω R the initial value Ω min q (t)/ζ(t) ≃ 29 at s = 2 is reduced to about 2.1 at s = 40 for γ = 0.1 and from about 12 at s = 2 to about 0.5 at s = 40 for γ = 0.95 (Ω min q (t)/ζ(t) ≃ 1 is reached at about s ≃ 29). Despite the improvement in the adiabaticity condition between the three different ramps, there are very little changes in Fig. 8, especially in atom number fluctuations. Reaching the limit where the ratio Ω min q (t)/ζ(t) ≫ 1 is both numerically and experimentally demanding in the case of optical lattices with large occupation numbers and a large number of sites. For instance, with the present parameter values, maintaining Ω min q (t)/ζ(t) 10 during the entire ramp to s = 40 for γ = 0.95 would already require a very slow ramp speed of δ ≃ 0.2 × 10 −3 ω R . The difficulty of adiabatically turning up a lattice potential in the experiments can severely limit possibilities to reach the superfluid Mott-insulator transition in the case of large occupation numbers [47]. Achieving a strong reduction in atom number fluctuations in the case of many atoms has been experimentally challenging even in the lattice systems of only a few sites [29].
The excitation of higher modes [ Fig. 7] induced during the ramping process indicate a nonadiabatic turning up of the lattice. The lattice-induced excitations are also reflected in relative atom populations between the two species (spin-1/2 waves). We show such excitations in the stable regime γ < 1 in Fig. 9 by displaying the overlap integral, defined by Eq. (25), between the wave functions of two atomic species for various values of inter-species nonlinear parameter γ = 0.1, 0.25, 0.5, 0.75. In deeper lattices the population difference is clearly increased. The effect of interactions, however, is again reduced as γ → 1.

Effects of temperature
Finite temperature in non-equilibrium quantum dynamics introduces additional noise in the system, increasing the atom number fluctuations. In experiments on atom number squeezing and reduced on-site atom number fluctuations in optical lattices the finite temperature has been an important factor limiting the achievable spin squeezing and the suppression of the atom number fluctuations [29,33]. Here we demonstrate the effects of the initial temperature of the atoms on the coherence properties of the twospecies system as the lattice potential is turned up. The temperature can be incorporated in the stochastic sampling of mode populations according to Eq. (32) in which case the width of the Gaussian stochastic distribution for the sampling of the initial state is increased due to thermal population of each phonon mode. In Fig. 10   on the relative intra-species phase coherence between the atoms in different lattice sites and on the on-site atom number fluctuations. The inter-species interaction parameter γ = 0.5. The corresponding stationary averaged values of atom number fluctuations n ( j) sqz and the relative intra-species phase coherence between the atoms in the adjacent lattice sites C 1 are shown in Fig. 11.

Unstable regime
In the previous section we considered a two-species condensate mixture in the regime where the spatial overlap of the two species is dynamically stable, corresponding to the values of the inter-species interaction strength γ 1. When the parameter γ is increased the system becomes dynamically unstable and the normal mode frequencies of Eq. (11) exhibit nonvanishing imaginary parts, indicating perturbations that grow exponentially in time. The instability criteria in different regimes for static and moving condensates were analyzed in detail in Refs. [8,16] and also the effects of matter-wave grating of the other species [60] have been studied. Phase separation is a generic phenomenon that occurs in different forms of matter. The phase separation instability condition for a two-species BEC in a lattice is analogous to the phase separation instability criterion of the two BEC components that occurs in free space when the square of the inter-species interaction coefficient exceeds the product  of the intra-species interaction coefficients.
Here we consider the unstable regime of (γ > 1) by first evaluating the thermal equilibrium state of the atoms in the initial state for some value of γ < 1 corresponding to a stable regime of overlapping two-species mixture. Stable initial ground-state configuration allows us to evaluate the statistical noise for the initial state of the TWA simulations within the Bogoliubov approximation. We then continuously turn up the optical lattice potential as in the dynamically stable case and immediately after the final lattice height is reached we change the inter-species interaction γ to the unstable regime. Varying the final lattice height and the atom number provides then information about the dependence of the phase separation dynamics on the lattice parameters and on quantum fluctuations. Experimentally the manipulation of the scattering lengths for BECs in order to drive the system from the stable to unstable regime has been realized for a 85 Rb-87 Rb condensate mixture using a Feshbach resonance [2]. Related experiments on single-component BECs by rapidly changing the scattering length from the stable positive to unstable negative value have generated a condensate collapse and the formation of bright solitons [61,62,63]. In two-species condensates the effective interaction strengths have been manipulated between stable and unstable regimes using dressed atomic states by electromagnetically-induced Raman transitions between the internal atomic states [35].

Domain formation
The system develops instability when the inter-species interaction exceeds a threshold γ ≃ 1 determined by the intra-species interactions. After the interaction parameters are switched to the unstable regime at the end of the lattice ramping the atom densities of the two BECs show violent phase separation dynamics and individual sites become dominantly occupied by single species alone. In the ground-state configurations of phase-separated systems one species typically forms a shell around the other one, minimizing the surface area between the two components [38]. In an optical lattice system we consider, quantum noise and the nonlinear interactions spontaneously break the symmetry of the uniform spatial configurations of the two species. Individual stochastic realizations of the TWA dynamics that represent possible outcomes of individual experimental runs show density domain formations of the two species. The system settles down to a metastable configuration of several interleaved spatially-separated density domains of the two species where the entire system may consist of multiple domain boundaries along the lattice. Individual sites are typically dominantly occupied by one atomic species alone, except close to the domain boundaries for the case of large domain length. In Fig. 12 we show typical single realization results for the atom density distribution for lattice heights s = 5, 10, 15, 20 after the system has phase-separated. The number of domains increases as the lattice height increases. The observed phenomenon is related to the experiments on a harmonically trapped 85 Rb-87 Rb condensate mixture in the absence of a lattice [2] in which case one of the species was found to split into multiple separated atom cloudlets that appeared as distinct holes in the density distribution of the other species. The condensation experiment was theoretically analyzed in Ref. [38] and it was demonstrated how the dynamics leads to continuous separation of the two species into smaller and smaller domains. As shown in Ref. [16], the formation of only a few density spikes and holes in the phase separation dynamics can be identified as a spontaneous generation of bound pairs of dark and bright solitons [12]. The phase separation provides a mechanism for the background of a dark soliton (density hole) in one species to stabilize a bright soliton (density spike) in the other species due to an effective trap that results from the repulsion between the two species [13]. In the present system the stabilization of the small domains is similar to the stabilization in energetically metastable particle-like solitons [13] and results in metastable configurations that are not necessarily energetically close to the ground state. The formation of the metastable states is a non-equilibrium process and the interleaved pattern does not represent a thermal state. One should compare the observed phase-separated state to the ground state that is a maximally phase-separated state and minimizes the surface area between the two components by forcing them to the opposite sides of the trap (one component to the right and the other one to the left).
In the spontaneous pattern formation, due to an instantaneous switch of the interaction strengths, the domain length is expected to be approximately determined by the wavelength of the phonon mode with the largest imaginary part of the eigenfrequency, since this corresponds to the unstable eigenmode that grows most rapidly (provided that the perturbation is strong enough to populate this mode). We can calculate the wavenumber of the mode q max analytically from the expressions of Ω q . The value of q max represents the fastest growing mode and gives an order-of-magnitude estimate for the domain length by l ∼ 1/|q max |. For J 1 = J 2 we obtain for the phase separation domain length where ∆ 2 12 > ∆ 11 ∆ 22 [where ∆ i j is defined in Eq. (8)] and the latter equality is valid for ∆ 11 = ∆ 22 . The domain length from Eq. (46) depends on the lattice height that modifies the dispersion relation. In deep lattices and with stronger interactions the domain length becomes shorter according to Eq. (46), as also numerically demonstrated in Ref. [39]. Note that reaching the ground state of maximally phase-separated state would require a very slow transition to the unstable regime so that only the lowest energy unstable mode is seeded in the process.
In Fig. 13 we show the classical field-theory estimates for the domain lengths obtained from Eq. (46) with l ∼ 1/|q max | for different values of the nonlinearity as a function of the final lattice depth. The domain length rapidly decreases as the lattice becomes deeper, eventually saturating due to the finite range of available q values in the lowest energy band.
In order to study the effect of quantum fluctuations of the atoms on the domain formation we vary the strength of quantum fluctuations in the simulations. We may continuously interpolate from the regime of strong quantum fluctuations to the classical mean-field limit by keeping the nonlinear interaction strengths χ i j constant, but varying the atom number [64]. This is done by changing the parameter Ω ⊥ a j j /(Nd) for constant χ i j . In the Bogoliubov approximation the nonlinearities χ i j fix the number of excited-state atoms depleted from the ground state, so varying the total atom number changes the depleted fraction due to quantum  20 to 1000, with L = 128. The limit n → ∞ corresponds to the classical mean-field limit. We find that both the number of domains and the fluctuations in the number of domains are increased due to quantum fluctuations. The domain boundaries between the two species may be viewed as defects and by introducing an electromagnetic coupling between the two components that mixes the atom populations provides a phase-separation scheme suitable for testing the Kibble-Zurek mechanism for defect formation in phase transitions [42,43]. Changing the value of the interaction parameter through the phase transition point spontaneously breaks the symmetry of the system in which case the formation of the defects is expected to depend on the timescale of the transition, providing an interesting link to condensed matter physics and cosmology.

Relative phase coherence and number fluctuations
In the unstable regime there is a dramatic loss in the relative phase coherence between the atoms in different lattice sites and the growth of on-site atom number fluctuations. The enhanced on-site atom number fluctuations may be understood in terms of the domain formation where the spatial location of the domain boundaries fluctuates from one stochastic realization to another. In Fig. 15  we show the atom number fluctuations in one lattice site for one of the atomic species (left column) and the relative intra-species phase coherence between the atoms in adjacent sites (right column) for different lattice heights s = 10, 20, 30, 40 and different inter-species interactions γ = 1.5 (top row) and γ = 2.0 (bottom row). Even though there is a rapid loss of phase coherence due to the dynamical instability, the nonvanishing value of C 1 indicates that even the dynamically unstable system exhibits nonvanishing phase coherence. The loss of phase coherence is much faster due to dynamical instability than due to the ramping of the lattice. The stationary averaged values as a function of the lattice height after the system has reached a metastable configuration are shown in Fig. 16. We also display in Fig. 16 the density overlap integral between the two BEC components τ(t) [Eq. (26)] for different values of the inter-species interaction strength γ that indicates the degree of spatial phase separation between the components. In Fig. 17 we show the populations in the lowest five momentum modes in the unstable regime. The decay of population in the ground state is fast due to the dynamical instability. The rapid decay of the population is a characteristic feature of the instability. We also show the populations of all the momentum states at a given time (t = 20000/ω R and γ = 1.1) for a single stochastic realization and for an ensemble average over many realizations.

Concluding remarks
We demonstrated via numerical simulations how a simple approximate stochastic phasespace method can provide valuable information about a two-species superfluid system in the presence of significant quantum fluctuations. We identified contributions of quantum fluctuations in the pattern formation of the two-species system that results from the modulational instability of dynamically unstable excitations. The parameter space of the twocomponent condensate system in the lattice is especially large. Although the essential effects of the phase separation dynamics, the loss of the relative phase coherence between the atoms and the reduced atom number fluctuations were captured by the selected parameter regimes, the parameter space could be explored in more detail. In particular, novel phenomena would be observed in the case of moving condensates. The centre-of-mass motion of the atoms results in dissipative transport properties [65,64]. The stability criteria are also changed when the velocities of the two BECs are different and one of the BECs is not in the normal dispersion regime [8].
Moving condensates may be experimentally studied in a combined optical lattice and a harmonic trap by suddenly displacing the harmonic trap, e.g., by using a magnetic field gradient, in which case the displacement excites dipolar oscillations of atoms along the lattice direction [65]. The other alternatives to create an analogous effect, e.g., are to use a movingstanding wave, so that the atoms experience a moving optical lattice potential [66] introduce a phase shift for the hopping amplitudes of the atoms between adjacent sites, or to make the hopping amplitudes time-dependent by a periodically pulsating lattice.
It would be particularly intriguing to investigate soliton structures in a two-species BEC. Modifying the ratio χ 22 /χ 11 [16], e.g., by unbalanced relative populations, or inducing flow instabilities [17,67] would lead to the emergence of vector solitons. The solitons can be persistent (metastable) or pulsating. Another potential application follows from the observation that repulsive interactions and quantum fluctuations lead to suppressed atom number fluctuations in the lattice. Two-species systems have already been used in experimental realizations of spin-squeezed atom interferometry [24] and squeezing can be employed in interferometers to reach sub-shot-noise accuracies that are not achievable using classical interferometers limited by the standard quantum limit [68,69,70,71]. n ( j) gr ≃ N j /L denotes the ground state atom number per site of the species j and w ( j) q,η ≡ u ( j) q,η −v ( j) q,η . We may introduce the corresponding phase operator at the site l aŝ q,η e iql − H.c. , (A.3) for which the commutator [n ( j) l ,φ ( j) l ] = i and we have defined r ( j) q,η ≡ u ( j) q,η + v ( j) q,η . Then the on-site atom number fluctuations in the lth site (∆n ( j) l ) 2 and the relative phase fluctuations between the atoms in the kth and lth site (∆ϕ ( j) kl ) 2 , respectively, read 1 N j q,η=± |r ( j) q,η | 2 sin 2 q(k − l) 2 (2n q,η + 1) = 1 N j q 1 2ǫ q (2n q,− + 1)Ω − q + (2n q,+ + 1)Ω + q sin 2 q(k − l) 2 , (A. 5) wheren q,± is the thermal population of the phonon mode (q, ±) in the lattice given in Eq. (32).
At T = 0 we haven q,η = 0 and we can evaluate the expressions analytically by replacing in the continuum limit the momentum sums by integrals For the atom number fluctuations we obtain In the nonlinear limit of (∆ 11 ± ∆ 12 ) ≫ J this simplifies to Similarly, for the relative phase fluctuations between the atoms in the adjacent sites (k − l = 1) we obtain In the nonlinear case (∆ 11 ± ∆ 12 ) ≫ J this reads (∆ϕ ( j) l,l+1 ) 2 = 1 n ( j) gr π √ 2J ∆ 11 + ∆ 12 + ∆ 11 − ∆ 12 . (A.10) Similar relationships as Eq. (A.9) may also be calculated for different values of k in (∆ϕ ( j) l,k ) 2 in the continuum limit.
The analytic expression (A.9) for the relative phase fluctuations can be used to evaluate as displayed in Fig. A1 in which case we show the relative phase coherence between the atoms in the adjacent sites C ′ 1 together with the on-site atom number fluctuations. We also use the continuum limit approximation to calculate (φ k+ j −φ k ) 2 for different values of j in order to obtain the coherence along the lattice C ′ j as a function of the site separation j. The parameters of Fig. A1 are the same as those used in in the TWA simulations for the dynamically stable T = 0 cases, with the dependence of the hopping amplitude J and the nonlinearity χ 11 = χ 22 on the lattice height determined by Eqs. (36) and (37), respectively, where we set the value χ 11 = 0.6ω R at s = 2. The atom number N = 2560 and the number of sites L = 64.
In Fig. A1 the relative phase coherence between the atoms in the adjacent sites decreases rapidly as a function of the lattice depth and increases as the inter-species interaction strength γ is increased closer to the onset of the phase-separation instability. The numerical values of the nearest-neighbour coherence are very close to those of the TWA simulations in Fig. 1, but the long-range coherence values are higher than in the TWA case [ Fig. 2]. One should note, however, that C ′ 1 does not include the atom number contributions incorporated in the definition of C 1 [Eq. (41)] which is used in analyzing the relative phase coherence in the TWA numerics.
We can implement a nonlinear least square fit for the coherence along the lattice C ′ j using a trial function