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

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 non-equilibrium 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.

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 ⁸⁷Rb in a ring of circumference $L=96$ $\mu$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 ⁸⁷Rb, but that $a_{12}$ is approximately half the naturally occurring value, making our system strongly immiscible with $\Delta=0.25$. We will discuss this issue further in Sec. 9. Furthermore, we take the total number of atoms to be $N=5\cdot 10^{4}$ and the transverse trapping frequency $\omega_{\perp}=2\pi\cdot 2$ kHz which sets the spin healing length $\xi_{\rm s}=\hbar/\sqrt{2m\rho g_{\rm s}}$ to be 1.5 times the transversal size of the system — enough to freeze the transversal spin degrees of freedom. The spin interaction constant is given by $g_{\rm s}=(g_{11}+g_{22}-2g_{12})/2$. We simulate phase transitions with quench times $\tau_{\rm 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)=|\psi_{1}(x)|^{2}-|\psi_{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_{\rm d}$ versus the quench time $\tau_{\rm Q}$ is plotted in Fig. 2(b). We fit the power law $N_{\rm d}=a_{\rm uni}/\tau_{\rm Q}^{n}$ to the results of the simulations with $\tau_{\rm Q}>2$ ms for which we simulated 1000 trajectories to minimise statistical uncertainty. We find a scaling exponent of $n=0.341\pm 0.006$ while $a_{\rm uni}=8.63\pm 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 $\nu=1/2$ and $z=1$ as derived in A and in Ref. [49].

The simulation data in Fig. 2(b) shows a deviation from the expected scaling behaviour for fast quench times $\tau_{\rm Q}<2$ ms. As quenches become faster ($\tau_{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 $\xi_{\rm 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_{\rm d}^{\rm max}\approx L/\xi_{\rm 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_{\rm d}$ for these quenches in Fig. 2(d). We note that for $\tau_{\rm 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_{\rm 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.

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\omega_{\rm x}^{2}x^{2}/2$ with $\omega_{\rm x}=2\pi\cdot 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 $|\psi_{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)=|\psi_{1}(x)|^{2}-|\psi_{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 $\gamma$ that limits the counting region for domains with a particle density larger than the threshold value $\rho_{\rm cut}=\gamma\rho(0)$. Figure 3(b) shows the behaviour of the average number of formed domains in function of the quench time for several values of $\gamma$. The scaling exponent is again determined by fitting a power law to the data points for quench times where we simulated 1000 trajectories. Figure 3(c) plots the extracted scaling exponent $n$ as a function of $\gamma$. Here we can see that $n$ shows a dependence on the counting region with values ranging from $n=0.289\pm 0.019$ to $n=0.473\pm 0.012$, but is a relatively constant $n\approx 0.47$ for $\gamma>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 $\gamma>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 $\tau_{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_{\rm d}$.

The emergence of a different scaling exponent for $\gamma>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 $\Omega_{\rm cr}(x)$ and consequently of the control parameter $\epsilon_{\rm 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.

While the simulations of the experiment in the harmonic trap suggest a scaling exponent for the number of defects of $n\sim 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 $\rho(x)$ is the Thomas-Fermi density of the trapped BEC system we are “simulating” [Eq. (15)], and the circumference of the ring $L=110$ $\mu$m $<2R_{\rm TF}$ (equivalent to $\gamma=0.15$) so that we avoid any divergence of $\Omega(x,t)$ at the edge of the simulated system where $\rho(x)\rightarrow 0$. In this situation the propagating phase transition front is due to the spatially dependent coupling parameter $\Omega(x,t)$, rather than the spatially dependent density. By design, this simulation has the same control parameter for the quench as $\epsilon_{\rm 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\pm 0.015$. Thus we conclude that the exponent of $n\sim 0.5$ that was determined in the harmonic trap is robust.

The success of the quantum “simulation” in Sec. 5 suggests that we can potentially engineer a spatially-dependent coupling $\Omega(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 $\Omega(x,t)=\max[0,2\Omega_{\rm cr}^{0}(x)(1-t/\tau_{\rm Q})]$, where $\Omega_{\rm cr}^{0}(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 $\Omega(x,t)$ to be proportional to the critical coupling at any time

where we have used the dependence of $\Omega_{\rm cr}(x)$ on the density. We have also introduced the function $c(t)=2(1-t/\tau_{\rm Q})$ which describes the proportionality between $\Omega(x,t)$ and $\Omega_{\rm cr}(x)$. Eqs. (28–29) constitute a system of equations describing the mutual relation between the condensate density $\rho(x,t)$ and the coupling $\Omega(x,t)$ that we set to a time-varying multiple of the critical coupling $\Omega_{\rm cr}(x,t)$ decided by $\rho(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 $\int dx\rho(x)=N/2$ we obtain an equation for $\Omega(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$ $\mu$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 $\tau_{\rm 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 $\tau_{\rm 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\int_{\rm-R_{\rm TF}}^{\rm R_{\rm TF}}\rho(x)dx$, we obtain the chemical potential

and the Thomas-Fermi radius

From (32) it can be seen that $R_{\rm TF}$ increases as $c(t)\to 1^{+}$, explaining the breathing in the case of a coupling quench from $c(0)>1$ to $c(t_{\rm final})=0$.

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_{\rm F}$. We control the velocity of the front by using the following quench of the coupling

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_{\rm F}$ is shown in Fig. 9, which indicates a value for the critical speed of $v_{\rm F}/v_{\rm s}\approx 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 $\xi_{\rm s}$ hence it is unlikely to have substantial effects.

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 time $\hat{t}$ 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 $\tau_{Q}$ that is different on each side of the critical point. We choose a quench

For a given quench time $\tau_{\rm Q}$, from Eq. (34) we can see that the rate of change of $\Omega$ 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 $\tau_{\rm Q}$ for this asymmetric quench. We fit the data with a power law $N_{\rm d}=a_{\rm asy}/\tau_{\rm Q}^{n}$ and derive the factor $a_{\rm asy}=8.74\pm 0.06$ which is in good agreement with the same quantity extracted from the data of Sec. 3 in Fig. 2(c) of $a_{\rm uni}=8.63\pm 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\hat{\xi}$, 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_{\rm asy}\approx 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.

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 $\Delta=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 ⁸⁷Rb and ²³Na naturally satisfy the criteria of naturally strong immiscibility, but the combination of the $|F=1,m_{\rm F}=+1\rangle$ and $|F=2,m_{\rm F}=-1\rangle$ hyperfine states in ⁸⁷Rb has an interspecies Feshbach resonance [68] that can be used to tune $\Delta$ to $0.8$ while keeping $g_{11}\approx g_{22}$. Although the use of a Feshbach resonance enhances losses, smaller values of $\Delta$ accelerate the process of domain formation allowing for the observation of the density pattern. We estimate that for a ring BEC with $L=50$ $\mu$m and $N=5000$, $\Delta=0.8$ yields $\xi_{\rm S}=0.4$ $\mu$m, $\tau_{0}=351$ ms and $N_{\rm d}^{\rm max}\approx 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 ⁸⁷Rb and was able to measure the phase separation of the dressed states across the critical point by changing the coupling strength $\Omega$. 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 ($\Delta\ll 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.

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 $\Omega$ 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.