Causality and defect formation in the dynamics of an engineered quantum phase transition in a coupled binary Bose-Einstein condensate

Continuous phase transitions occur in a wide range of physical systems, and provide a context for the study of non-equilibrium dynamics and the formation of topological defects. The Kibble-Zurek (KZ) mechanism predicts the scaling of the resulting density of defects as a function of the quench rate through a critical point, and this can provide an estimate of the critical exponents of a phase transition. In this work we extend our previous study of the miscible-immiscible phase transition of a binary Bose-Einstein condensate (BEC) composed of two hyperfine states in which the spin dynamics are confined to one dimension [J. Sabbatini et al., Phys. Rev. Lett. 107, 230402 (2011)]. The transition is engineered by controlling a Hamiltonian quench of the coupling amplitude of the two hyperfine states, and results in the formation of a random pattern of spatial domains. Using the numerical truncated Wigner phase space method, we show that in a ring BEC the number of domains formed in the phase transitions scales as predicted by the KZ theory. We also consider the same experiment performed with a harmonically trapped BEC, and investigate how the density inhomogeneity modifies the dynamics of the phase transition and the KZ scaling law for the number of domains. We then make use of the symmetry between inhomogeneous phase transitions in anisotropic systems, and an inhomogeneous quench in a homogeneous system, to engineer coupling quenches that allow us to quantify several aspects of inhomogeneous phase transitions. In particular, we quantify the effect of causality in the propagation of the phase transition front on the resulting formation of domain walls, and find indications that the density of defects is determined during the impulse to adiabatic transition after the crossing of the critical point.


Introduction
The study of non-equilibrium quantum systems is exemplified by systems undergoing a quantum phase transition, in which the ground state undergoes a macroscopic change of symmetry as a Hamiltonian parameter is changed through a critical value [1]. These phase transitions often result in the formation of topological defects, occurring when the topology of the degenerate ground states of the new phase is non-trivial [2]. In this case, particular classes of excited states of the system, the topological defects, cannot be continuously transformed to the ground state and their existence is therefore protected by topological conservation laws. We can potentially unveil some of the mysteries behind out-of-equilibrium systems by studying the mechanism of formation of topological defects. These are imprinted in the system during the "turbulent" era of a phase transition but survive its subsequent evolution. Like fossils, topological defects provide information about a period that cannot be easily modelled analytically. The ability to engineer a quantum phase transition and recover information from topological defects is an invaluable resource for the study of nonequilibrium systems [3,4].
The Kibble-Zurek mechanism is a theory developed by Kibble [5] and Zurek [6] with the goal of predicting the scaling of the density of topological defects formed with the rate at which the critical point of a phase transition is crossed. The theory is based on the assumption that every system undergoing a phase transition experiences a phase of impulsive behaviour where the topological defects are "crystallised" in the system. The breakdown of adiabaticity, due to the divergence of the response time near a critical point, is a characteristic of every continuous phase transition and the theory offers a model for the formation of defects across a wide range of systems, from cosmology [7] to condensed matter [8,9,10,11,12,13,14,15,16]. While initially applied to classical phase transitions occurring at finite temperature [5,6], more recently the scenario has been formulated and successfully applied to quantum phase transitions at zero temperature [17,18,19,20,21,22,23,24].
The development of isolated, yet highly tunable ultra-cold gas experiments have provided exciting opportunities for the study of quantum matter and quantum dynamics. Recent advances in the production and manipulation of Bose-Einstein condensates (BEC) offer the possibility of observing a variety of quantum phase transitions, from the study of the Ising model with ultra-cold gases [25] to the exploration of the Mott insulator to superfluid phase transition [26]. Experimentalists in degenerate quantum gases have access to a large set of tools like Feshbach resonances [27], arbitrary trapping [28] and single-atom imaging [29] that offer them the ability to engineer Hamiltonians and to control their parameters. Furthermore the production of the first optically trapped BEC [30] led to the exploration of the richer physics of multi-component condensates. Our ability to accurately manipulate and prepare a wide range of states of condensates with internal degrees of freedom allowed the observation of spin changing dynamics [31] and gave us access to the phase diagram of spinor BECs [32]. Multi-component BECs were also the key piece for the first observation of spin-orbit coupling in neutral atoms [33].
In this work we extend a proposal for an experiment to test the KZ theory using a binary BEC consisting of a two hyperfine states of a single atomic species [34]. Binary condensates can be classified as miscible or immiscible, depending on whether the two states can naturally coexist in space [35]. However, it has been shown that introducing a coupling between the two states can transform an immiscible condensate into a miscible one, and that there is a quantum phase transition between the two states [36,37]. In the strong coupling regime the mean-field ground state of the system consists of a superposition of the two states whose density then resembles the miscible phase. Here we consider a naturally immiscible condensate beginning in the ground state of the strong coupling regime, which is subsequently driven back to the immiscible phase by quenching the coupling to zero (see Fig. 1) [34]. If the quench is non-adiabatic, regions of the system that are not causally connected will undergo the phase transition independently and the final configuration is a random pattern of domains of the two components, Fig. 1(c). The KZ mechanism is demonstrated by counting the number of domains formed as a function of the characteristic quench time. When applied to a BEC in a ring trap [28,38,39] this scheme is an experimentally feasible candidate to provide a definitive demonstration of the KZ mechanism [34]. In our study we perform our simulations using the truncated Wigner approximation (TWA) [40] to numerically simulate the quantum dynamics of the BEC during the quench. This phase-space method maps quantum operators to stochastic variables [40] to include quantum corrections to the mean-field dynamics, allowing us to capture the symmetry breaking that occurs in the phase transition.
In the vast majority of theoretical and experimental works, the KZ mechanism has been studied for homogeneous phase transitions [6,17,14], i.e. phase transitions in which the critical point is traversed simultaneously for all spatial regions of the system. However, ultra-cold quantum gases are often confined to trapping potentials which lead to spatially varying phases and inhomogeneous phase transitions [41]. In such situations additional considerations, such as the effects of causality, can change the quantitative details of the phase transition, and the KZ theory needs to be modified [42,43,44]. In the implementation we describe we can study the effects of inhomogeneities by engineering the phase transition in uniform systems. Spatially shaping the coupling quench, for example, allows the possibility of "simulating" an inhomogeneous phase transition in systems with otherwise straightforward dynamics such as a homogeneous 1D ring BEC. The precise level of control required over the phase transition can be experimentally achieved with the use of modern experimental techniques that, for example, offer the possibility of shaping the spatial intensity profile of a laser using holography [45,46]. Engineered Hamiltonian quenches are potentially a powerful tool to study the characteristics of phase transitions in a controlled fashion that can be applied to theoretical and experimental studies of other quantum phase transitions. This paper is structured as follows. In Sec. 2 we introduce the model of binary condensates and outline the principles behind the KZ mechanism. In Sec. 3 we model the miscible-immiscible phase transition in a ring BEC, and demonstrate that the number of defects scales with the quench time as predicted by the KZ theory [34]. We then consider the same phase transition in a harmonically trapped BEC in Sec. 4. We find that the counting of domains can be problematic, and we find a different scaling exponent with respect to the uniform case. In Sec. 5 we engineer a spatial quench of the coupling capable of "simulating" the phase transition of an harmonically trapped BEC in a ring geometry. This engineered phase transition allows us to verify the exponent observed in the previous section. In Sec. 6 we extend the method of spatially shaping the coupling to invert the process, i.e. we attempt to simulate a homogeneous phase transition in an harmonically trapped system. We quantify the effects of causality for an inhomogeneous phase transition by accurately identifying the subsonic and supersonic regimes of the phase transition in Sec. 7. In Sec. 8 we demonstrate that for this system the domain walls are seeded after the passing of the critical point, and finally in Sec. 9, we explore the experimental feasibility of our scheme and summarise our results. In the strong coupling regime state (a) becomes miscible. Quantum noise has been added that results in the formation of domains (see text). (c) Quenching the coupling Ω(t) to zero brings the system back to its immiscible phase. If the quench is nonadiabatic a random pattern of domains is created and the system is left in an excited state. We propose to test the Kibble-Zurek mechanism by counting the number of formed domains in function of the quench time.

Binary Bose-Einstein condensates
Several authors have described the physics of binary BEC mixtures [35,47] and coupled binary BECs [48,49]. We define the Hamiltonian for our coupled 1D system aŝ whereĤ 0 ,Ĥ I , andĤ C are the single-particle, interaction, and coupling Hamiltonians respectively, defined aŝ The Bose field operatorψ i (x) for component i annihilates a particle at position x, and obeys the commutation relations The transversal dimensions, assumed to be harmonically trapped with frequency ω ⊥ , have been integrated out from the three-dimensional Hamiltonian, yielding the 1D interaction constants g ij = 2h 2 a ij /(ma ⊥ ), where a ij is the scattering length and a ⊥ = h/mω ⊥ . The coupling HamiltonianĤ C depends on the coupling strength Ω(t) and on the detuning δ between the light-field and the energy difference of the states. Here we assume the coupling is on resonance and set δ = 0 for the remainder of this paper.
The nature of the ground state of a binary system is determined by the balance between the interaction constants g ij [35]. If the intraspecies interactions dominate, g ii g 12 , the energy is minimised by lowering the individual densities of both components, which then occupy all the available volume. In this miscible phase, the two species coexist at every point. However if the interspecies interaction dominates g 12 g ii , the energy of the system is minimised by phase separation, which minimises the contribution of the g 12 n 1 (x)n 2 (x) term. We refer to this phase as the immiscible phase. By defining then we have [35] ∆ > 1 miscible phase, However, the introduction of the coupling Ω between the states can qualitatively change this picture. In the strong coupling regime, when Ω is much larger than a threshold value Ω cr , the ground state of the system approaches the superposition state |ψ = (|1 + |2 )/ √ 2. In this situation the atoms of a naturally immiscible mixture then spatially coexist as for the miscible phase. The calculation of the value of the critical coupling Ω cr can be found in Appendix A.

Kibble-Zurek mechanism
To formulate the KZ mechanism [5,6] for this system, we define the dimensionless control parameter which quantifies the distance of the system from the critical point Ω cr . For a continuous phase transition in the thermodynamic limit, both the equilibrium correlation length ξ and relaxation time τ of the system diverge as where ν and z are the critical exponents of the transition and define its universality class. The constants ξ 0 and τ 0 depend on the particular characteristics of the system, such as the atomic species, particle density, trapping frequencies, and so on. The correlation length represents the average size of the regions where the system has uniform characteristics, while the relaxation time characterises the time needed by the system to adjust to an external change. For a system dynamically approaching a critical point, there exists a moment when the relaxation time is equal to the characteristic time scale on which the environment is changing. The KZ mechanism predicts that beyond this time the system is not able to follow the quench and remain in quasi-equilibrium -instead it undergoes impulsive behaviour which freezes the configuration of defects in space. We consider a quench of the coupling parameter for our system of the form where τ Q is the quench time. The moment when the relaxation time τ [ (t)] becomes larger than the characteristic time of the quench | (t)/˙ (t)| determines the freezing timet through the condition Solving Eq. (10) with the aid of Eqs. (7)-(8) for the freezing time giveŝ The correlation length at freezing time,ξ = ξ 0 (τ Q /τ 0 ) ν/(1+νz) , fragments the system in a series of regions with uniform characteristics. Since topological defects can form only on the boundaries between these regions [50] the number of formed defects is given by

Simulation method
Our goal in this paper is to perform a computational study of the dynamics of the engineered quantum phase transition. In order to do so, we must go beyond the mean-field approximation for the system. For a system with a uniform initial density, dynamically crossing the critical coupling Ω cr leads to a modulational instability -such that the initial state is dynamically unstable. In mean-field simulations such instabilities are seeded by numerical noise. An approximate method for simulating beyond mean-field methods for ultra-cold gases is the truncated Wigner approximation [51,52,40]. Briefly, this is a phase space method that simulates the evolution of the Wigner function for a system by means of stochastic trajectories. It is approximate in that the stochastic representation neglects some terms in the equation of motion involving third and higher orders derivatives of the Wigner function with respect to the phase space variables. However, these have a small contribution for short times when the number of particles in the system is much larger than the number of modes that are simulated [52,40].
The net result is that the equations of motion for the system are simply the appropriate Gross-Pitaevskii equations for the coupled condensates, but the initial state must be sampled from the Wigner distribution. For a weakly interacting BEC at zero temperature, this corresponds to populating the Bogolioubov quasiparticle modes with half a particle of classical noise, which represents the initial quantum fluctuations. Expectation values of symmetrised products of equal time quantum field operators are calculated by the averages of the corresponding phase space variable over a number of trajectories beginning with different initial conditions. In our simulations, the initially small "quantum noise" then seeds the modulational instability, and this leads to macroscopically different outcomes as might be expected to be realised in an experiment. Indeed, it has been plausibly argued that individual trajectories can be roughly interpreted as being equivalent to single experimental runs [40]. This is the approach that we take for the analysis in this paper.
The equations of motion for the two components We construct the initial states for our simulation for the homogeneous system beginning from the even superposition ψ where N is the total number of particles and L is the length of the ring. For the case of an elongated BEC the initial Gross-Pitaevskii solution is found by imaginary time propagation of Eq. (13) in the strong coupling regime. As prescribed by the TWA, the initial state is then formed by generating the complex Gaussian noises η i with variance η * i η j = δ ij /2, and forming ψ( where ] are the amplitudes of the i-th Bogoliubov mode, found by solving the Bogoliubov-de Gennes equation for a binary BEC [37,53].

Homogeneous Bose-Einstein condensate in a ring trap
We begin our investigation of the KZ mechanism in the engineered phase transition in a binary BEC by considering a system confined in a ring trap, as was originally described in Ref. [34]. For our simulation we choose to simulate 87 Rb in a ring of circumference L = 96 µm, and take the scattering lengths to be a 11 = a 22 = a 12 /2 = 1.325 nm. We note that a 11 and a 22 are close to the measured values for 87 Rb, but that a 12 is approximately half the naturally occurring value, making our system strongly immiscible with ∆ = 0.25. We will discuss this issue further in Sec. 9. Furthermore, we take the total number of atoms to be N = 5 · 10 4 and the transverse trapping frequency ω ⊥ = 2π · 2 kHz which sets the spin healing length ξ s =h/ √ 2mρg s to be 1.5 times the transversal size of the systemenough to freeze the transversal spin degrees of freedom. The spin interaction constant is given by g s = (g 11 + g 22 − 2g 12 )/2. We simulate phase transitions with quench times τ Q over three order of magnitudes, in the range [0.1, 125] ms. The large ratio between the number of particles N and the number of simulated modes M = 1024 ensures the validity of the truncated Wigner approximation [52].
The time evolution of the condensate density is shown in Fig. 2(a) where it is possible to observe an example of the final pattern of domains. The number of domains formed at the end of the simulation is identified as the number of zero crossings of the function f (x) = |ψ 1 (x)| 2 − |ψ 2 (x)| 2 . This method of domain counting can easily be adopted experimentally by performing absorption imaging of the two components after Stern-Gerlach separation [54,55]. The mean number of formed domains N d versus the quench time τ Q is plotted in Fig. 2(b). We fit the power law N d = a uni /τ n Q to the results of the simulations with τ Q > 2 ms for which we simulated 1000 trajectories to minimise statistical uncertainty. We find a scaling exponent of n = 0.341 ± 0.006 while a uni = 8.63 ± 0.03. The scaling exponent n is in good agreement with the theoretical prediction of the KZ theory n = 1/3, obtained from Eq. (12) with the mean-field critical exponents ν = 1/2 and z = 1 as derived in Appendix A and in Ref. [49].
The simulation data in Fig. 2(b) shows a deviation from the expected scaling behaviour for fast quench times τ Q < 2 ms. As quenches become faster (τ Q decreases), the KZ mechanism predicts a decreasing correlation length at the freezing time, and hence a smaller average domain size. In contrast, in our scheme the correlation length has a lower bound given by the spin healing length ξ s , inducing the saturation of the number of domains for fast quenches as seen in Fig. 2(b). The value of the saturation we extract from our data is compatible with the predicted maximum number given by N max d ≈ L/ξ s . However, we also observe a mean number of domains significantly higher than the prediction for fast quenches [dashed box in Fig. 2(b)]. We plot the time evolution of N d for these quenches in Fig. 2(d). We note that for τ Q < 1.5 ms the number of domains is still decreasing at the end of the integration time. It is in principle possible to extend the integration time to extract the true value of N d after stabilisation occurs and study phase ordering [71] but we encountered significant numerical error and the "thermalisation" of the quantum noise when pushing the integration time beyond 600 ms [52]. The latter is a known signature of the invalidity of the TWA for long evolution times, and is particularly accentuated by fast coupling quenches. These result in the rapid growth of the population of short wavelength modes, leaving the system in a highly excited state that enhances the thermalisation process.

Inhomogeneous Bose-Einstein condensate in a elongated harmonic trap
We now move to consider the same experiment performed for a binary BEC in an elongated harmonic trap, as is common to many experimental setups for ultra-cold gases. We take V (x) = mω 2 x x 2 /2 with ω x = 2π · 5 Hz, with all other parameters the same as in Sec. 3. An example trajectory for the trapped case is shown in Fig. 3(a).
For the homogenous BEC it was relatively straightforward to identify the location of domain walls. For the inhomogeneous system, near the boundary of the condensate the densities |ψ i (x)| 2 are of a comparable magnitude to the noise that is added to the initial state, and the density difference function f (x) = |ψ 1 (x)| 2 − |ψ 2 (x)| 2 can develop zeros that are not clearly domain walls. Experimentally an equivalent issue arises due to the presence of thermal and instrumental noise in the imaging system as discussed in [56] for quantum vortices.
To attempt to address this issue, we introduce a new parameter γ that limits the counting region for domains with a particle density larger than the threshold value ρ cut = γρ(0).  to n = 0.473 ± 0.012, but is a relatively constant n ≈ 0.47 for γ > 0.3, pointing to a systematic difference of the scaling exponent for an harmonically trapped BEC with respect to the homogeneous case. For large counting regions γ > 0.25 (where domains are counted in the low density region of the condensate), it can be seen in Fig. 3(b) the number of domains for larger quench times τ Q seems to converge to a value larger than predicted by a scaling law. This effect can be traced to the miscounting of noise as domains which results in an overestimate of the physical number of domains N d . The emergence of a different scaling exponent for γ > 0.3 in the harmonically trapped case can be linked to the effects of inhomogeneity and causality, and this is the major topic we wish to address in this paper. Qualitatively, the value of the critical coupling is connected to the particle density, therefore the phase transition in an harmonically trapped BEC becomes inhomogeneous. The centre of the BEC, the point of highest density, enters the new phase first and the transition proceeds then towards the edges, as can be seen in Fig. 3(a) where the domain formation near the edges occurs at later times, and a propagating phase front can clearly be identified. In a trapped system the spatial dependence of Ω cr (x) and consequently of the control parameter tr (x, t) breaks the translational symmetry giving a preferred direction of motion for the newly formed domains. In systems of finite size this has the effect of increasing the annihilation rate of domains or the rate at which they escape our counting regions. This effect becomes increasingly important for slower quenches, and increases the observed scaling exponent.

Causality in the harmonically trapped system
The issue of causality is related to the speed of the front separating regions in the unbroken (old) and broken (new) symmetry phase. It has been proposed [42,57,58] that in the case where the front is slower than the speed of sound, information can propagate from regions in the new phase to regions in the old phase. This exchange of information influences the choice of symmetry of the system in the old phase, reducing the number of defects that are eventually formed. To analyze this phenomenon, we start by deriving an expression for the density of a binary condensate in the strong coupling regime. Applying the Thomas-Fermi (TF) approximation for the density [59] we have for x smaller than the Thomas-Fermi radius R TF defined by the solution of µ − V (R TF ) + hΩ(0) = 0, where µ is the chemical potential resulting from N = 2 dxρ(x). As described above, the critical coupling Ω cr is now spatially inhomogeneous, with dependence given bȳ hΩ cr (x) = (g 12 − g)ρ(x). We can modify our definition of the control parameter Eq. (6) to reflect the spatial inhomogeneity of the phase transition where we have used Ω(t) = max[0, 2Ω cr (0)(1−t/τ Q )]. In inhomogeneous phase transitions, different regions of the system experience different effective quench time as a consequence of the spatially dependent critical coupling. The effective quench time is again defined as the change rate of the newly defined control parameter tr (x, t), We note that the effective quench time τ Q (x) approaches zero for x → R TF . Therefore, outer regions of the BEC are subject to a smaller correlation length at freezing time, and the average size of the domains decrease in the proximity of the condensate edges, see Figs. 3(a) and 4. Solving the equality tr (x F , t F ) = 0 for x F defines the trajectory of the front which in the case of harmonic potential moves according to We have verified that Eq. (18) accurately reproduces the trajectory of the front by comparing it with the simulations, as shown in Fig. 4. The speed of the front during the transition follows from Eq. (18) as

Kibble-Zurek scaling exponent for the harmonically trapped system
In order for the new phase to transfer information to the old phase, the front speed has to be smaller than the relevant speed of sound at freezing time given byv =ξ/τ = 1+νz . For the phase transition under consideration here (z = 1), v is independent of the quenching time and equal to the local speed of sound v s (x) = (g + g 12 )ρ(x)/(2m) [60,37]. Solving the inequality for x defines the regionX within which the formation of defects is suppressed. For the parameters used in this work Eq. (20) has no solution for τ Q < τ cr Q ≈ 150 ms. Thus for a large range of quench times we simulate, the front is always faster than the sound, and the phase transition is not affected by the broken symmetry choices in the neighbouring regions. For τ Q > τ cr Q on the other hand, the speed of the front equals the speed of sound in two points (see Fig. 5a, where the pink solid line compares the speed of the phase transition front to the local speed of sounds for τ Q = 223 ms.) When the front first enters the condensate at the centre it moves faster than the sound for a short distance, until it slows into a subsonic regime. Approaching the edges, where the speed of sound goes to zero, the front again becomes supersonic. As the size of the region with suppressed domain formation increases with the increasing quench time, the effect of causality on inhomogeneous phase transitions is to increase the scaling exponent, by further reducing the total number of domains formed compared to the fully supersonic case. We can compute the number of formed domains in this scenario by integrating the effective quench time τ Q (x) over the supersonic region: wherex 1 andx 2 are respectively the inner and outer boundaries of the subsonic regionX. When integrated numerically with ν = 1/2 and z = 1, Eq. (21) yields a scaling exponent n ≈ 0.37.
In more detail, we note that the comparison of the local speed of sound to the front propagation speed does not provide a rigorous criterion for delimiting the suppression of domain formation. In fact, the authors of Ref. [57,58] note that defect formation can occur for v F /v s ≥ 0.5. However, the introduction of any proportionality constant in Eq. (20) modifies the value of the critical quench time τ cr Q , whereas the value of the scaling exponent remains n ≈ 0.37. We address the issue of the critical value for the speed of the front in Sec. 7.

Spatial density of defects
We turn now to the analysis of the dependence of the spatial density of defects on the condensate density for a harmonically trapped BEC. A comparison between the front speed for a range of quench times and the local speed of sound is shown in Fig. 5(a). In Fig. 5(b) we plot the mean spatial distribution of the defect density at the end of the integration time for the same quench times as in Fig. 5(a). We can see that for fast quenches (for which the front speed is always larger than the speed of sound) the domain density increases towards the edge of the condensate as suggested by Eq. (17). The rise in defect density for x > 0.8R TF is due to the difficulty of distinguishing noise from domains in the low density region, as described in Sec. 4. Figure 5(b) shows that for slower quenches, the density of domains begins to decrease for increasing x, rather than increasing as is observed for faster quenches. This is likely due to the suppression of domain formation as described above. The onset of the suppression is expected to result in a sudden and sharp decrease in the number of domains. However, we have been unable to observe this clearly in the simulation results, as the positions of the domains are not fixed when Ω(t) > 0, and so domains initially seeded in the supersonic region are able to move into the subsonic region before they become clearly distinguishable.

Inhomogeneous phase transition simulated in the ring geomety
While the simulations of the experiment in the harmonic trap suggest a scaling exponent for the number of defects of n ∼ 1/2, this is still open to question. To confirm the modified scaling exponent derived in Sec. 4 we now "simulate" the inhomogeneous phase transition of the harmonically trapped BEC. We acheive this by performing numerical simulations of a binary BEC in a ring trap subject to a spatially dependent quench of the coupling strength designed to reproduce the physics of the trapped system. We use where ρ(x) is the Thomas-Fermi density of the trapped BEC system we are "simulating" [Eq. (15)], and the circumference of the ring L = 110 µm < 2R TF (equivalent to γ = 0.15) so that we avoid any divergence of Ω(x, t) at the edge of the simulated system where ρ(x) → 0.
In this situation the propagating phase transition front is due to the spatially dependent coupling parameter Ω(x, t), rather than the spatially dependent density. By design, this simulation has the same control parameter for the quench as tr (x, t) in Eq. (16). However, we avoid the domain-counting problem described in Sec. 4 as the BEC has constant density around the ring. The evolution of the density of one component of the BEC is shown in Fig. 6(a). The critical coupling is first reached at x = 0, and the front of the phase transition moves around the ring in exactly the same manner as in the trapped BEC with a spatially constant coupling. The scaling of defect number for this simulation is shown in Fig. 6(b), and we find a scaling exponent of n = 0.497 ± 0.015. Thus we conclude that the exponent of n ∼ 0.5 that was determined in the harmonic trap is robust.

Kibble-Zurek scaling exponent for the inhomogeneous phase transition
We now estimate the Kibble-Zurek scaling exponent expected for the inhomogeneous phase transition in a ring BEC. For a ring BEC the condition Eq. (20) for domain suppression due to causality becomes where on the right hand side we have inserted the constant speed of sound in the ring [37]. The number of domains scales as function of the quench time as wherex is the solution to the equality of Eq. (23). If the size of the system L is such that the inequality Eq. (23) is never satisfied, the temporal dependence in Eq. (24) factors out. This will increase the number of domains compared to the case of the homogeneous transition, but leaves the scaling exponent unchanged at n = 1/3. However, if there is a region where the speed of the phase transition front is less than the speed of sound, the defect-suppression region grows with the quench time and the scaling exponent is n = 4/3.
For the results we present in Fig. 6(b) the quench times are small enough that there is no region of domain suppression. However, we find the scaling exponent to be n ≈ 0.5, in disagreement with the prediction above. We believe this discrepancy is caused by the breaking of translational invariance, as we describe below in Sec. 5.2.

Translational invariance and domain annihilation
For a BEC in a ring trap the system is symmetric under rotations about the center of the ring. This is reflected in the fact that the Hamiltonian Eq. (1) of the main text is invariant under translations x → x + α when V (x) = 0 (the definition of a homogeneous system). However, harmonically trapped condensates and spatially inhomogeneous quenches of the coupling break this symmetry. To clearly demonstrate the effect of breaking translational invariance on the phase transition we simulate both a homogeneous and an inhomogeneous quench in a uniform system, but seed the domain formation "by hand." We begin with the initial state where k is an integer that determines the number of seeded domains, and A = 0.1 is the amplitude of the seed. The coupling is for the homogeneous quench and for the inhomogeneous quench emulating the phase transition in the harmonically trapped BEC. Both choices ensure that Ω(x, 0) < Ω cr everywhere at t = 0. The domains will grow from the seed with no propagating front for the phase-transition, but the nonzero coupling will affect the domain dynamics. We plot typical results in Fig. 7. For fast quenches we find that the domains do not have time to move or decay, and the number of domains observed at the end of the integration time is the same as the number seeded for both homogeneous and inhomogeneous quenches. However, for slow quenches we find that Ω I (x, t) results in a systematically smaller number of domains than Ω H (x, t) for the same quench time, due to the motion and subsequent annihilation of domain walls. This effect will increase the scaling exponent for the number of domains, and for our parameters it increases it to n ≈ 0.5.

Simulating a Homogeneous Phase Transition in an Harmonically Trapped BEC
The success of the quantum "simulation" in Sec. 5 suggests that we can potentially engineer a spatially-dependent coupling Ω(x, t) such that the phase transition in an harmonically trapped BEC happens everywhere simultaneously -that is, the phase transition is homogeneous. We describe our efforts to realise this below.
A harmonically trapped BEC reacts to a coupling quench by changing its shape. This means that a coupling quench of the form Ω( where Ω 0 cr (x) is the initial critical coupling, will not produce an homogeneous quench as it fails to take into account the change of shape of the system density. We extend Eq. (15) to take into account a spatially and time dependent coupling with the following ansatz for the density This reflects Eq. (15) in the case where the coupling is spatially inhomogeneous and uses the local density approximation, i.e. the condensate density is influenced only by the local value of the trapping potential and coupling between the components. We also assume the coupling Ω(x, t) to be proportional to the critical coupling at any timē where we have used the dependence of Ω cr (x) on the density. We have also introduced the function c(t) = 2(1 − t/τ Q ) which describes the proportionality between Ω(x, t) and Ω cr (x).
Eqs. (28-29) constitute a system of equations describing the mutual relation between the condensate density ρ(x, t) and the coupling Ω(x, t) that we set to a time-varying multiple of the critical coupling Ω cr (x, t) decided by ρ(x, t). The condition expressed by Eq. (29) ensures the simultaneous crossing of the critical point for c(t) = 1. Substituting Eq. (29) into Eq. (28) and solving the equation dxρ(x) = N/2 we obtain an equation for Ω(x, t) for an harmonically trapped system where G(t) = g 12 + g − c(t)(g 12 − g). We have implemented the coupling quench of Eq. (30), and the evolution of the density of one of the components is shown in Fig. 8. We find that the phase transition is effectively homogeneous in the range |x| < 50 µm. However, in deriving Eq. (30) we have assumed that the condensate adiabatically follows the spatial dynamics of the quench at any time, i.e. the quench time is much larger than the condensate reaction time. For faster quenches the condensate fails to follow the changes in the coupling and the conditions leading to Eq. (30) break down. As τ Q decreases, the region where the transition is effectively homogeneous rapidly shrinks or even disappears, complicating the process of domain counting. This leaves us with a narrow range of quench times τ Q with a homogeneous transition, which is insufficient to reliably extract a scaling exponent. From Fig. 8, it is also evident that the quench of Eq. (30) results in the onset of breathing dynamics that further complicates the counting process. In order to explain the origin of the breathing mode we derive the Thomas-Fermi radius of the condensate in the approximation given by Eqs. (28)(29). Inserting the solution of the system of equations into the normalisation condition N = 2 RTF −RTF ρ(x)dx, we obtain the chemical potential and the Thomas-Fermi radius From (32) it can be seen that R TF increases as c(t) → 1 + , explaining the breathing in the case of a coupling quench from c(0) > 1 to c(t final ) = 0.

Threshold Value for Causality
Previous analysis of the effect of causality on the formation of defects found that defect suppression occurs when the speed of the front is a fraction of the speed of sound [57,58].
In this section we determine the value of the critical velocity of the phase front below which the formation of domains is suppressed for the coupled binary BEC. To do so we simulate the inhomogeneous phase transition in a ring with a front moving with the velocity v F . We control the velocity of the front by using the following quench of the couplinḡ where the initial distance between the fronts w is introduced to avoid the effect of the front entering the system. We simulate this system using the same parameters as Sec. 3 with w = L/20 for a time of t = 150 ms. The mean number of defects formed as a function of the speed of the front v F is shown in Fig. 9, which indicates a value for the critical speed of v F /v s ≈ 0.37. In the mean field approximation the transition between the supersonic and subsonic regime is expected to occur at a precise fraction of the speed of sound. However the presence of fluctuations in the density of the system means that there is some uncertainty in the the local speed of sound, smearing out the exact value of the transition. In Ref. [58] the authors determine the relevant speed for the onset of causality using a moving step function for the control parameter. In this work we use a spatially continuous control parameter and consequently the phase transition occurs on a finite width. This finite width is however small and comparable with the spin healing ξ s hence it is unlikely to have substantial effects.

When is the Defect Density Determined?
Historically the KZ mechanism has been at the centre of a debate about the timing of the birth of defects [61,62]. On one hand, it has been suggested that the appearance of defects can be traced to the fluctuations that freeze out at a timet before the transition occurs, i.e. when the dynamics first changes from being adiabatic to impulsive [63]. However, there is also a symmetric point after the transition, when the system dynamics again becomes adiabatic. An alternative viewpoint has been expressed that it is this point, after the transition has occurred, when the eventual defect density is determined [62,70]. In our numerical simulations we can probe the time at which the domain density is determined by employing a piecewise linear coupling quench for the ring-shaped system, with a quench time τ Q that is different on each side of the critical point. We choose a quench For a given quench time τ Q , from Eq. (34) we can see that the rate of change of Ω before the critical point is 3/2 times smaller compared to the critical point. The choice of the relatively small factor of 3/2 is due to the finite range of quenching times we can reliably simulate. Very fast quenches result in a larger number of domains that relaxes even after the end of the quench; on the other hand, very slow quenches are more likely to experience problems with the truncated Wigner approximation (Sec. 3). We performed simulations of this system with the parameters otherwise as described in Sec. 3, and in Fig. 10 we plot the number of domains formed as function of τ Q for this asymmetric quench. We fit the data with a power law N d = a asy /τ n Q and derive the factor a asy = 8.74 ± 0.06 which is in good agreement with the same quantity extracted from the data of Sec. 3 in Fig. 2(c) of a uni = 8.63 ± 0.03.
The KZ mechanism does not provide a conclusive prediction for the total number of topological defects formed in a phase transition, instead focusing on their scaling. This uncertainty is usually treated by saying that the typical distance between defects is fξ, where f is unknown and typically varies from 1-10 [64]. However, in simulating the quenches in Fig. 10, we have only altered the slope of the ramp before the critical point, while all the other parameters remained identical. As a result it seems reasonable to assume that the pre-factor f is the same for both scenarios.
With the assumption that f is fixed, if the defects were decided before the critical point, the absolute number of domains for the scenario with quench time 3/2 times slower should have been reduced by a factor of (3/2) n = 1.14, where we have used n = 1/3. Then we would have expected to measure a prefactor of a asy ≈ 7.54. This is well outside the uncertainty range of the value that we have measured. The equivalence of the total number of domains for both type of quenches strongly suggests the defect density is decided after the transition.
This conclusion is compatible with the physical mechanism of domain formation in our system. In the miscible-immiscible quantum phase transition, domains form as a consequence of modulational instability of the system [65]. In this situation the population of unstable modes with complex eigenvalues grows exponentially after the transition, suggesting that the number of formed domains is decided by the characteristic length of the most unstable mode at the moment of the return of adiabaticity.

Experimental Feasibility
We now examine the experimental feasibility of our scheme. The miscible and immiscible phases were experimentally observed in various binary systems [66] and recently both phases were realised with the same pair of atomic species using Feshbach resonances [67]. The results we present in this paper are obtained in the regime where the two components are strongly immiscible with ∆ = 0.25. This means that the system spin healing length is relatively short, and leads to both a large number of domains and their straightforward identification. While such a strongly immiscible system may be difficult to realise in practice, it does not present an overwhelming obstacle to the experimental realisation of the scheme described here.
No pair of hyperfine states of 87 Rb and 23 Na naturally satisfy the criteria of naturally strong immiscibility, but the combination of the |F = 1, m F = +1 and |F = 2, m F = −1 hyperfine states in 87 Rb has an interspecies Feshbach resonance [68] that can be used to tune ∆ to 0.8 while keeping g 11 ≈ g 22 . Although the use of a Feshbach resonance enhances losses, smaller values of ∆ accelerate the process of domain formation allowing for the observation of the density pattern. We estimate that for a ring BEC with L = 50 µm and N = 5000, ∆ = 0.8 yields ξ S = 0.4 µm, τ 0 = 351 ms and N max d ≈ 50.
Recently, Lin et al. reported the observation of the miscible-immiscible quantum phase transition in spin-orbit coupled Bose-Einstein condensates [33]. The group coupled two Zeeman sub-levels of the F = 1 hyperfine state of 87 Rb and was able to measure the phase separation of the dressed states across the critical point by changing the coupling strength Ω. The method is not affected by increased inelastic losses, as outlined in the previous paragraph for Feshbach resonances, and is capable of achieving higher separation regimes (∆ 1) that make the task of counting domains easier. Similar results have also been reported by Nicklas et al. [69]. However, we note that the spatial configuration of the dressed state is not directly accessible, but is instead reconstructed from absorption imaging of the bare components. Hence, the higher separation and stability come at the price of a more complicated detection process for the number of domains. In addition, the results presented in Ref. [69] and our own numerical simulations suggest that the appearance of domains in the dressed state is slower than in the procedure proposed here. As a consequence, this alternative scheme could potentially result in higher domain loss before the counting, further altering the scaling exponent.

Conclusions
In conclusion, we have shown that the number of defects formed in the coupling induced miscible-immiscible phase transition in a ring-shaped binary BEC scales as predicted by the Kibble-Zurek theory. The scaling exponent we determine numerically for the number of defects for this system agrees with that predicted by using the mean-field critical exponents. The results suggest the scheme is a good candidate for the experimental testing of the predictions of the KZ mechanism in an experiment with ultra-cold gases, taking advantage of the system's isolation from the environment and the precise control of the parameter Ω that can be provided by microwave or laser coupling.
We have also examined the phase transition in a binary BEC in an elongated harmonic trap, and found that this yields a modified scaling exponent compared to the homogeneous case. We have shown how an engineered coupling quench can be used to simulate inhomogeneous phase transitions in a ring BEC, and derived the scaling exponent of the harmonically trapped case independent of the domain counting problem inherent to trapped BEC systems.
Engineered quenches can provide a systematic way to study the relatively unexplored problem of inhomogeneous phase transitions. Using this technique we have shown how it is possible to invert the process and realise a homogeneous phase transition over a large spatial region of an harmonically trapped BEC. Although we were unable to drive the entire condensate through the critical point at the same time due to the breakdown of our approximation, the approach clearly establishes a path to simplify phase transitions in inhomogeneous systems.
Finally we have addressed two specific issues of the KZ mechanism: the relation between causality and defect formation, and the critical moment at which the defect density is decided. For the former we have derived a threshold value for the velocity of a front at which defect formation is suppressed in our system. For the latter, using a temporally asymmetric quench, we have supplied additional evidence in support of the case made in Ref. [62], i.e. that the density of defects is decided after the transition for this particular realisation of the KZ mechanism.
where ξ B =h/ √ 2mgρ, and ρ = N/L is the linear density. The diagonalisation of the matrix S results in a pair of eigenvalues Γ i = 2{1+1/ √ ∆, 1−1/ √ ∆+hΩ/gρ} and a transformation matrix U such that A = U SU −1 where A is a diagonal matrix. In the basis (φ 1 , φ 2 ) T = U (φ 1 , φ 2 ) T , Eq. (A.3) takes the form We can then extract information about the miscible-immiscible phase transition by examining the effective correlation lengths ξ i = ξ 2 B /Γ i . As expected, an uncoupled mixture (Ω = 0) has a critical point for ∆ = 1. We are interested here in the value of the critical coupling which turns an immiscible condensate into a miscible one. From ξ 2 it follows that when ∆ < 1 the critical point is It is also possible to quantify the critical exponent ν by noticing that the effective correlation length scales as ξ 2 = ξ B /Γ 2 = ξ B /| | 1/2 , which combined with Eq. (7) implies ν = 1/2. Similar results for the critical coupling and critical exponent were obtained in Ref. [49] from the energy spectrum of the system. We quote here only the energy gap E gap ∝ h Ω(t)[Ω(t) − Ω cr ] of the spectrum in the mean-field approximation. The presence of a gapped spectrum indicates that the phase transition is continuous. The relaxation time is related to the energy gap through the relation τ =h/E gap ∝ 1/| | 1/2 which implies the dynamical critical exponent is z = 1.