Defect production in nonequilibrium phase transitions: Experimental investigation of the Kibble-Zurek mechanism in a two-qubit quantum simulator

Systems passing through quantum critical points at finite rates have a finite probability of undergoing transitions between different eigenstates of the instantaneous Hamiltonian. This mechanism was proposed by Kibble as the underlying mechanism for the formation of topological defects in the early universe and by Zurek for condensed matter systems. Here, we use a system of nuclear spins as an experimental quantum simulator undergoing a non-equilibrium quantum phase transition. The experimental data confirm the validity of the Kibble-Zurek mechanism of defect formation.


I. INTRODUCTION
When a system is driven through a continuous phase transition at a finite rate, domain structures can arise due to spontaneous symmetry breaking. With the growth of the domains, they can approach each other and generate, e.g., topological defects. This process was initially proposed by Kibble for studying the cosmological phase transition in the early universe [1,2], and developed by Zurek in condensed matter systems [3][4][5]. Today, it is known as the Kibble-Zurek mechanism (KZM) and has become a universal theory for studying non-equilibrium dynamics in both classical and quantum systems [6]. The KZM owes its appeal and broad applicability to its universal and scale invariant power law prediction, based on simple arguments about second-order phase transitions. Typical examples include critical phenomena in the quantum Ising model [7,8]. The KZM has been supported by experiments in physical systems, such as superfluids [9,10], superconductors [11,12], Bose-Einstein condensates [13,14], colloidal monolayers [15], ion crystals [16,17], and more solid materials [18,19].
The interest in the KZM was recently enhanced by its relation with the established Landau-Zener (LZ) model [20], where the critical point can be modeled with a sim- * Electronic address: jingfu@e3.physik.tu-dortmund.de † Electronic address: fernando@cucchietti.com ‡ Electronic address: laflamme@iqc.ca § Electronic address: Dieter.Suter@tu-dortmund.de ple avoided level crossing in a two level system. In a true second-order phase transition, which can only occur in many-body systems, we would observe instead a symmetry-breaking accompanied by a closing of the gap. The LZ model cannot be as complex as a large quantum many-body system with a phase transition, although we expect to capture at least qualitatively the dynamics of the change of properties during a quantum phase transition (QPT). Similar to a second-order QPT, the defect density created during the passage through the critical point can be controlled through the quench rate versus the inverse of the gap; if this is large, it results in a high defect density. In the opposite limit of a slow quench, the transition through the critical point becomes almost adiabsatic and the defect density tends to zero. The LZ model provides a good test-bed for studying the KZM in a well-controlled system, such as by quantum simulations [21][22][23][24].
The essential concept behind the KZM can be introduced as follows. A thermodynamic system initially in thermal equilibrium evolves under some control parameter, e.g., pressure or temperature, that drives a phase transition. If the control parameter changes slowly enough compared with the relaxation time τ of the system, the evolution is adiabatic, i.e., the system can adjust to the new conditions so that it remains in thermal equilibrium. However, if the control parameter changes too rapidly, the system cannot adjust to the new conditions sufficiently fast, and it is driven out of equilibrium. This process happens mostly near the critical point of the phase transition, where the relaxation time τ diverges. In this regime, the dynamics of the system becomes impul-sive and its state is effectively frozen, i.e., the change of the control parameter in the evolution only introduces an overall factor to the waveform function of the system. Following the KZM model [3], we approximate the evolution as discontinuous between the two regimes, and define the boundary as the freeze-out timet, where the relaxation time of the system is equal to the time required to reach the critical point at the given (constant) scan rate.
Existing experimental investigations of the KZM by quantum simulators [21][22][23][24] relied on single-qubit systems. In the case of the Ising model, it is possible to map the dynamics to the LZ model [22][23][24]. In this article, we present an experimental quantum simulation of the KZM in a system of two interacting qubits, where the defect generation occurs non-locally, simultaneously creating quantum entanglement. In our experiment, we exploited two spin-qubits in a nuclear magnetic resonance (NMR) system, using the natural Ising interaction of the spins. A small transverse field was applied to generate the level-anticrossing for simulating the second order QPT driven by a control field along z-axis. The defect generation and the freezing-out time were clearly demonstrated. Our implementation in a two-qubit system can be generalized to larger qubit systems, providing a testbed for simulating the KZM in many-body systems.

II. THEORETICAL MODEL
We first consider the phase diagram of an Ising model consisting of qubits 1 and 2. Its Hamiltonian is where σ i z denotes the z components of the Pauli operators and B z the magnetic field component and the coupling strength between the two qubits has been set to unity. By solving the Hamiltonian (1), we obtain the energy levels as 1+2B z , −1, and 1−2B z for the triplet eigenstates |00 , |φ + = (|01 + |10 )/ √ 2, and |11 , respectively. Here, we do not include the singlet state |φ − = (|01 − |10 )/ √ 2, since the symmetry of the Hamiltonian confines the evolution of a triplet state to the triplet subsystem. Depending on the control parameter B z , the ground state of the system is then .
|φ + is one of the Bell states, a maximally entangled state and therefore a useful resource in the field of quantum information [25,26]. As a function of the control parameter B z , the system undergoes a QPT with critical points B z = ±1 [27,28].
To simulate the KZM, we drive transitions of the system through the critical points at B z = ±1. These transitions can be adiabatic only if there is no exact crossing. These avoided crossings can be generated by adding a small magnetic field along the x-axis. The Hamiltonian of the system is then Figure 1(a) illustrates the anti-crossing energy levels for B x = 0.1. The instantaneous ground state can be represented as Analytical expressions for the coefficients c 0 , c + , and c 1 are given, e.g., in [27]. We initialize the system into its ground state at t = 0 and let the field B z change linearly with time, where k denotes the scan rate. At time t, the system has evolved into written in the eigenstates of the instantaneous Hamiltonian. The populations |a i | 2 depend on the rate at which the critical points are traversed. The propagator that converts |ψ g (0) into |ψ(t) can be represented as where T denotes the time-ordering operator. The dynamics of this system can also be simulated numerically. For these simulations, as well as for experimental simulations, we replaced the continuous scan by a stepwise constant effective Hamiltonian by dividing the total evolution period into N segments with duration δ. In the limit N → ∞ this becomes the ideal case. The evolution propagator is then where t = jδ with j = 0, 1..N . The initial state should be the ground state for B 0 → −∞, i.e. the state |00 . In the numerical simulations, we Figure 1(c) shows the populations |a i | 2 during the scan for B x = 0.1 and scan rates k = 1 and 1/20. At the beginning of the scan, when the gap between the ground and first excited states is large, the system can adjust to the change of the control parameter, and the evolution is therefore adiabatic. However, when the control parameter approaches the critical point, where the eigenstates change rapidly, the evolution becomes nonadiabatic [27,29] and the excited states become populated.
The characteristic time, at which the system can adjust to the change of the eigenstates, is the relaxation time τ . It can be defined by the inverse of the gap (in frequency units) between the ground and first excited states. Figure  1(b) shows the dependence of τ on the control field as the system goes through the phase transition region. It reaches a maximum at the two critical points in which the ground state is affected.
The non-equilibrium dynamics of this system can be put into the context of quantum critical systems [20]. Since the two critical points at B c z = ±1 are fully equivalent, we concentrate on the critical point at B c z = −1, i.e., B z < 0. If the system is initially in the ground state, the defect density D corresponds to the total population of the excited states,

III. EXPERIMENTAL PROTOCOL
As an experimentally accessible system for testing these predictions, we chose an NMR quantum simulator. As the quantum register, we used 13 C-labelled chloroform (CHCl 3 ) dissolved in d6-acetone. The proton and carbon nuclei are assigned as qubits 1 and 2. Data were taken with a Bruker DRX 700 MHz spectrometer. In the doubly rotating frame, the Hamiltonian of the NMR system is where ν denotes the offset, and J = 215 Hz denotes the coupling between the two qubits. The pseudo-pure initial state |00 was prepared by spatial averaging [30].
To optimize the experimental implementation, we used results from numerical simulations to choose experimental parameters and design the experiment protocol. The simulations showed: Here the discretization results in a reduction of the fidelity from 1 to 3. The defect density D remains close to zero (D < 0.005) for B z < −1.5, i.e. before the system approaches the critical point. We therefore started the scan at B z = −1.5. Figure 2 shows the pulse sequence for the experimental implementation. We divided the whole scan period into 15 segments with identical durations. The corresponding values of the control parameter B z are B z (jδ) = −1.5 + 0.1j. We first prepared the initial state P (0)|00 . For given fields B x and B z (t), the operator P (t) = e iβ(σ 1 generates the ground state |ψ g (B x , B z (t) from state |00 , where  Figure 2 shows the pulse sequence used for generating P (0), U (t), and P (t) † . The system was then allowed to evolve into the state U (t)P (0)|00 and the transformation P (t) † was applied, which converts the ground state of the final Hamiltonian into the |00 state. Then we applied a gradient pulse to eliminate coherence in this state since we only need diagonal terms. To obtain the overlap F (t), we performed partial quantum state tomography [33], using one π/2 read out pulse on qubit 1 and invoking the permutation symmetry of the two qubits.

IV. EXPERIMENTAL RESULTS AND DISCUSSION
We performed the experiment for different values of the transverse field, B x = 0.1 and 0.2 and scan rates k = 1, 1/2, 1/3 and 1/4. Figure 3 shows the results as circles.
The experimental results are complemented by two different types of simulations: the simulation of the ideal model is shown as a set of curves, while a simulation of the discretized experiment, using the actual experimental parameters but ideal pulses is shown as dots. The differ- ences between the simulated experiment and the theoretical curves reflect the approximations made by generating the stepwise constant propagator [see Eq. (12)]. Experimental and simulated data agree within the experimental uncertainties. Earlier benchmarking experiments [31] showed that the pulse control error is negligibly small. The differences between theory and experiment appear to be dominated by transverse relaxation (T 2 ) effects.
To estimate the effect of relaxation on the data, we performed numerical simulations for the case of the slowest scan rate (k = 1/4), using the experimentally determined relaxation times. The resulting defect densities are shown as blue squares in Figure 3. For a simplified description of the anticrossing region, we reduce the system to an effective two-level system described by the Hamiltonian Here, we have shifted the origin of the energy axis to the center of the gap. This model (16) is equivalent to the LZ model, which was used to simulate the KZM by Damski [20], up to a shift of B z (t). The quench time scale τ Q of the KZM is proportional to the ratio of the transverse field and the scan rate: τ Q = √ 2B x /k [20,34]. In Figure  4, we show the good agreement between the results by numerical simulation of the dynamics of H ef f and H in our region of interest.
Since we are mostly interested in the behavior near the critical point B c z = −1, we rescale the distance from the critical point by dividing the time |B z + 1|/k required to reach B c z by τ Q : ε = |B z (t) + 1|/(k τ Q ). In these units, the relaxation time becomes with the maximum of the relaxation time τ 0 = 1/(2 √ 2B x ). According to the KZM model, the boundary between the adiabatic and non-adiabatic ( or impulse [20]) regimes is given by the freeze-out timet, which is related to the  relaxation time τ by the equation [3,20] τ (t) = αt (18) where α = O(1) and independent of τ Q and τ 0 . The freeze-out timet indicates the distance to the critical point. The constant α is also related to the final defect density D f after the passage through the critical point The derivation of this relation is presented in the appendix.
Since we are only interested here in the effect of the passage through the first critical point, we determined the 'final' defect density D f at B z = −0.2, well after the first critical point, but before the second critical point. To test the scaling relation (19), we compare defect densities for a range of scan rates k and two transverse field strengths B x = 0.1 and 0.2. Figure 5 shows the dependence of D f on τ Q /τ 0 (or the scan duration scaled as B 0 /(4B 2 x ) ) for two sets of simulated data and for the experimental data of Figure 3. Linear fits of the simulated and experimental data sets yield the parameters listed in Table I. These results agree well with the theoretical predictions of Eq. (19) and α = O(1), as well as with the LZ-formula, which predicts α = π/2 [34,35]. The difference between simulated and actual experiments is mostly due to transverse relaxation.

V. ENTANGLEMENT FROM THE NON-EQUILIBRIUM PHASE TRANSITION
As a function of the control parameter, the two-qubit system used for these quantum simulations has different ground states; for |B z | > 1, the ground state is separable, for |B z | < 1 it is entangled, see Eq. (2). Accordingly, we prepared the system in a separable state, but depending on the adiabaticity, the passage through the critical points creates entanglement in the system. This entanglement generation can be controlled by the transverse field and the scan rate of the control field [36]. Figure 6 shows the results of numerical simulations, using the concurrence [37] as the entanglement measure. When B x is small and the scan is fast, such as in the case of the red dot-dashed curve in figure (a), the system remains in the separable state and almost no entanglement is generated, corresponding to many defects. For larger B x and slower  scan rate, the probability that the system remains in the ground state increases and the two qubits become entangled during the passage through the critical point, while the defect density decreases. In the limit of a slow scan, the evolution becomes adiabatic, and the entanglement becomes a measure of the QPT [27]. At intermediate scan rates, the concurrence shows significant oscillations, indicating that the system is far from equilibrium, in a superposition of the two lowest eigentstates. The lower sweep rates require longer scan times, e.g., about 30 and 90 ms for k = 1/10, and 1/30, respectively. In such cases, the errors from the relaxation effects increase. One possible solution for exploring this regime would be to choose a system with stronger couplings, such as dipolar coupled spins in liquid crystal solvents [38].

VI. CONCLUSION
We have experimentally simulated the KZM of defect generation in an interacting quantum system consisting of a pair of coupled nuclear spins. We controlled the spins by applying near-resonant radio-frequency fields. Our model contains several critical points; here, we focused on the one that is encountered first as the control parameter is swept from large negative to large positive fields. We demonstrated the formation of defects during the evolution starting from the ground state of the initial Hamiltonian. Like in the LZ model, the relaxation time of the system passes through a maximum at the critical point, but does not diverge, in contrast to the original KZM. We also used our model system to verify the validity of the Zurek equation and determine the freeze-out time.
Our work can be considered as a first step of a more general strategy to study the KZM in interacting many body systems with QPTs. While we used a comparatively simple two-qubit system, it can be readily extended to larger systems, which allow to simulate more complex quantum critical phenomena, such as in the quantum Ising model [39] or the generation of entanglement in the non-equilibrium phase transition.