Experimentally testing quantum critical dynamics beyond the Kibble–Zurek mechanism

The Kibble–Zurek mechanism (KZM) describes the dynamics across a phase transition leading to the formation of topological defects, such as vortices in superfluids and domain walls in spin systems. Here, we experimentally probe the distribution of kink pairs in a one-dimensional quantum Ising chain driven across the paramagnet-ferromagnet quantum phase transition, using a single trapped ion as a quantum simulator in momentum space. The number of kink pairs after the transition follows a Poisson binomial distribution, in which all cumulants scale with a universal power law as a function of the quench time in which the transition is crossed. We experimentally verified this scaling for the first cumulants and report deviations due to noise-induced dephasing of the trapped ion. Our results establish that the universal character of the critical dynamics can be extended beyond KZM, which accounts for the mean kink number, to characterize the full probability distribution of topological defects. Pioneered in the cosmological setting, the Kibble–Zurek mechanism (KZM) describes the universal dynamics across a phase transition, leading to the breakdown of adiabatic dynamics and the formation of topological defects. The authors present an experimental study of the universal critical dynamics of a 1D quantum Ising chain driven across the paramagnet-to-ferromagnet phase transition using a trapped-ion quantum simulator, and characterize the probability distribution of topological defects, which is found in excellent agreement with theoretical predictions.

T he understanding of nonequilibrium quantum matter stands out as a fascinating open problem at the frontiers of physics. Few theoretical tools account for the behavior away from equilibrium in terms of equilibrium properties. One such paradigm is the so-called Kibble-Zurek mechanism (KZM) that accounts for the universal critical dynamics of a system driven across a phase transition. Pioneering insights on the KZM were conceived in a cosmological setting 1 and applied to describe thermal continuous phase transitions [2][3][4] . The resulting KZM was later extended to quantum phase transitions [5][6][7][8] , see ref. 9 for a review. Its central prediction is that the average total number of topological defects 〈n T 〉, formed when a system is driven through a critical point in a time scale τ Q , is given by a universal power law hn T i $ τ Àβ Q . The power-law exponent β = Dν∕(1 + zν) is fixed by the dimensionality of the system D and a combination of the equilibrium correlation length and dynamic critical exponents, denoted by ν and z, respectively. Essentially, the KZM is a statement about the breakdown of the adiabatic dynamics across a critical point. As such, it provides useful heuristics for the preparation of ground-state phases of matter in the laboratory, e.g., in quantum simulation and adiabatic quantum computation 10 . It has spurred a wide variety of experimental efforts in superfluid helium [11][12][13] , liquid crystals 14,15 , convective fluids 16,17 , superconducting rings [18][19][20] , trapped ions [21][22][23] , colloids 24 , and ultracold atoms [25][26][27][28][29] , to name some relevant examples. This activity has advanced our understanding of critical dynamics, e.g., by extending the KZM to inhomogeneous systems 30,31 . In the quantum domain, experimental progress has been more limited and led by the use of quantum simulators in a variety of platforms [32][33][34] .
Beyond the focus of the KZM, the full counting statistics encoded in the probability distribution can be expected to shed further light. The number distribution of topological defects has become available in recent experiments 35 . In addition, the distribution of kinks has recently been explored theoretically in the onedimensional transverse-field quantum Ising model (TFQIM) 36 , a paradigmatic testbed of quantum phase transitions 37 . Indeed, it has been argued that the distribution of topological defects is generally determined by the scaling theory of phase transitions and thus exhibits signatures of universality beyond the KZM 36 .
We aim at validating this prediction experimentally using a single-trapped ion as a quantum simulator and probe the quantum critical dynamics of the TFQIM in momentum space. Specifically, we show that the measured distribution of kink pairs is Poisson binomial and that low-order cumulants scale as a universal power law with the quench time.

Results
Quantum critical dynamics. The TFQIM is described by the HamiltonianĤ that we shall consider with periodic boundary conditionsσ z 1 1 σ z Nþ1 and N even. Additionally,σ z m andσ x m are Pauli matrices at the site m and g plays the role of a (dimensionless) magnetic field that favors the alignment of the spins along the x-axis. This system exhibits a quantum phase transition between a paramagnetic phase (|g| ≫ 1) and a doubly degenerated phase with ferromagnetic order (|g| ≪ 1). There are two critical points at g c = ±1. We shall consider the nonadiabatic crossing of the critical point g c = −1 after initializing the system in the paramagnetic phase and ending in a ferromagnetic phase. Our quantum simulation approach relies on the equivalence of in the TFQIM in one spatial dimension and an ensemble of independent two-level systems. This result forms the basis of much of the progress on the study of quantum phase transitions and can be established via the Jordan-Wigner transformation 38 , Fourier transform, and Bogoliubov transformation 39 . We detail these steps in the Supplementary Note 1, where it is shown that the TFQIM Hamiltonian can be alternatively written in terms of independent modes as 37,40 whereγ k are quasiparticle operators, with k labeling each mode and taking values k ¼ π Other physical observables can as well be expressed in both the spin and momentum representations. We are interested in characterizing the number of kinks. As conservation of momentum dictates that kinks are formed in pairs, we focus on the kink-pair number operatorN The later can be equivalently written asN ¼ P k>0γ y kγk , wherê γ y kγk is a Fermion number operator, with eigenvalues {0, 1}. As different k-modes are decoupled, this representation paves the way to simulate the dynamics of the phase transition in the TFQIM in "momentum space": the dynamics of each mode can be simulated with an ion-trap qubit, in which the expectation value ofγ y kγk can be measured. To this end, we consider the quantum critical dynamics induced by a ramp of the magnetic field in a time scale τ Q that we shall refer to as the quench time. We further choose g(0) < −1 in the paramagnetic phase. In momentum space, driving the phase transition is equivalently described by an ensemble of Landau-Zener crossings. This observation proved key in establishing the validity of the KZM in the quantum domain 7,32,33 : the average number of kink pairs hN i ¼ hni after the quench scales as in agreement with the universal power law hni KZM / τ À ν 1þzν Q , with critical exponents ν = z = 1 and one spatial dimension (D = 1).
The full counting statistics of topological defects is encoded in the probability P(n) that a given number of kink pairs n is obtained as a measurement outcome of the observableN . To explore the implications of the scaling theory of phase transitions, we focus on the characterization of the probability distribution P (n) of the kink-pair number in the final nonequilibrium state upon completion of the crossing of the critical point induced by Eq. (3). Exploiting the equivalence between the spin and momentum representation, the dynamics in each mode leads to two possible outcomes, corresponding to the mode being found in the excited state (e) or the ground state (g), with probabilities p e = p k and p g = 1 − p k , respectively. Thus, one can associate with each mode k > 0 a discrete random variable, with excitation probability p k ¼ hγ y kγk i. The excitation probability of each mode is that of Bernoulli type. As the kink-pair number n is associated with the number of modes excited, P(n) is a Poisson binomial distribution, with the characteristic function 36 e PðθÞ ¼ Àππ associated with the sum of N/2 independent Bernoulli variables. The mean and variance of P(n) are thus set by the respective sums of the mean and variance of the N/2 Bernoulli distributions characterizing each mode, 〈n〉 = ∑ k>0 p k and Var(n) = ∑ k>0 p k (1 − p k ). The full distribution P(n) can be obtained via an inverse Fourier transform. The theoretical prediction of the kink distribution P(n) relies on knowledge of the excitation probabilities {p k }, which can be estimated using the Landau-Zener formula (details are given in Supplementary Note 2).
Experimental design. Experimentally, the excitation probabilities {p k } can be measured by probing the dynamics in each mode, that is simulated with the ion-trap qubit. The dynamics of a single mode k is described by a Landau-Zener crossing, that is implemented with a 171 Yb + ion confined in a Paul trap consisting of six needles placed on two perpendicular planes, as shown in Fig. 1a. The hyperfine clock transition in the ground state S 1∕2 manifold is chosen to realize the qubit, with energy levels denoted by 0 State preparation and manipulation. At zero static magnetic field, the splitting between 0 j i and 1 j i is 12.642812 GHz. We applied a static magnetic field of 4.66 G to define the quantization axis, which changes the 0 j i to 1 j i resonance frequency to 12.642819 GHz, and creates a 6.5 MHz Zeeman splittings for 2 S 1/2 , F = 1. In order to manipulate the hyperfine qubit with high control, coherent driving is implemented by a wave mixing method, see the scheme in Fig. 1c. First, an arbitrary wave generator (AWG) is programmed to generate signals~200 MHz. Then, the waveform is mixed with a 12.442819 GHz microwave (generated by Agilent, E8257D) by a frequency mixer. After the mixing process, there will be two waves~12.242 GHz and 12.642 GHz, so a high-pass filter is used to remove the 12.224 GHz wave. Finally, the wave~12.642 GHz is amplified to 2 W and used to irradiate the trapped ion with a horn antenna.
Measurement. For a typical experimental measurement in a single mode, Doppler cooling is first applied to cool down the ion 41 . The ion qubit is then initialized in the 0 j i state, by applying a resonant light at 369 nm to excite S 1/2 , F = 1 to P 1/2 , F = 1.
Subsequently, the programmed microwave is started to drive the ion qubit. Finally, the population of the bright state 1 j i is detected by fluorescence detection with another resonant light at 369 nm, exciting S 1/2 , F = 1 to P 1/2 , F = 0. Fluorescence of the ion is collected by an objective with a 0.4 numerical aperture. A 935 nm laser is used to prevent the state of the ion to jump to metastable states 41 . The initialization process can prepare the 0 j i state with fidelity >99.9%. The total error associated with the state preparation and measurement is measured as 0.5% (ref. 42 ).
In the Supplementary Note 1, we rewrite the explicit Hamiltonian of the k-th mode as a linear combination of a single two-levels system 7 . In this way, experimentally, the Hamiltonian of a single k-mode can be explicitly mapped to a qubit Hamiltonian describing a two-level system driven by a chirped microwave pulse, where Ω R = 4J/ℏ and Δ k ðtÞ ¼ 4J½gðtÞ þ cos k=ð_ sin kÞ are the Rabi frequency and the detuning of the chirped pulse, respectively. To measure the excitations probability after a finite ramp with normalized quench parameter A = Jτ Q /ℏ, we can vary g(t) from −5 to 0. For this purpose, a chirped pulse with length of T p ¼ 5AT R sin k and g(t) = 5t/T p − 5 is used in the experiment, with T R = 1/Ω R denoting the Rabi time. A typical operation process driven by the programmed microwave is illustrated on the Bloch sphere, see Methods section for details. First, the qubit is prepared in the ground state ofĤ TLS k at g(0) = −5. Then, g(t) is ramped from −5 to 0 with quenching rate 1/τ Q . Finally, the ground state of the finalĤ TLS k is rotated to qubit 0 j i, so that the excited state is mapped to 1 j i, which can be detected as a bright state. The measured excitation probability for different quench times is shown in Fig. 2.
Full probability distribution of topological defects. From the experimental data, the distribution P(n) of the number of kink pairs is obtained using the characteristic function Eq. (5) of a Poisson binomial distribution, in which the probability p k of each Bernoulli trial is set by the experimental value of the excitation probability upon completion of the corresponding Landau-Zener sweep. This result is compared with the theoretical prediction valid in the scaling limit-the regime of validity of the KZM-in which P(n) approaches the normal distribution, away from the adiabatic limit 36 A comparison between P(n) and Eq. (7) is shown in Fig. 3. The matching between theory and experiment is optimal in the scaling limit far away from the onset of the adiabatic dynamics or fast quenches, for which nonuniversal corrections are expected. We further note that the onset of adiabatic dynamics enhances non-normal features of the experimental P(n) that cannot be simply accounted for by a truncated Gaussian distribution, that takes into account the fact that the number of kinks n = 0, 1, 2, 3, … To explore universal features in P(n), we shall be concerned with the scaling of its cumulants κ q (q = 1, 2, 3… ) as a function of the quench time. We focus on the first three, theoretically derived in the scaling limit in the Supplementary Note 2. The first one is set by the KZM estimate for the mean number of kink pairs κ 1 (τ Q ) = 〈n〉 KZM . The second one equals the variance and as shown in the Supplementary Note 2 is set by Þhni KZM . Finally, the third cumulant is given by the third centered moment κ 3 ðτ Q Þ ¼ hðn À hniÞ 3 Indeed, all cumulants are predicted to be nonzero and proportional to 〈n〉 KZM (ref. 36 ), with the scaling of the first cumulant hni KZM / τ À 1 2 Q being dictated by the KZM. We compared this theoretical prediction with the cumulants of the experimental P(n) in Fig. 4, represented by dashed lines and symbols, respectively. As discussed in the Supplementary Note 2, deviations from the power law occur at very fast quench times, satisfying τ Q < ℏ/(π 3 J). However, this regime is not probed in the experiment, as all data points are taken for larger values of τ Q . Within the parameters explored, deviations at fast quenches with Jτ Q /ℏ~1 are due to the fact that the excitation probability p k in each mode is not accurately described by the Landau-Zener formula, as shown in Fig. 2. For moderate ramps, the power-law scaling of the cumulants is verified for q = 1, 2. The power-law scaling for The quantum critical dynamics of the one-dimensional transverse-field quantum Ising model is accessible in an ion-trap quantum simulator, which implements the corresponding Landau-Zener crossings governing the dynamics in each mode, labeled by wavevector k. The Rabi frequency Ω R = 4J/_ was set to 2π × 20 kHz in the experiment. For each value of the wavevector k, the excitation probability is estimated from 10,000 measurements. The shaded region describes the excitation probability predicted by the Landau-Zener formula. Deviations from the later become apparent for fast quench times, especially for large values of k. Error bars indicate the standard deviation over 10,000 measurements. Fig. 3 Probability distribution of the number of kink pairs P(n) generated as a function of the quench time. The experimental kink-pair number distribution for a transverse-field quantum Ising model with N = 100 spins is compared with the Gaussian approximation derived in the scaling limit, Eq. (7), ignoring high-order cumulants. The mean and width of the distribution are reduced as the quench time is increased. The experimental P(n) is always broader and shifted to higher kink-pair numbers than the theoretical prediction. Non-normal features of P(n) are enhanced near the sudden quench limit and at the onset of adiabatic dynamics. Error bars indicate the standard deviation over 10,000 measurements. Fig. 4 Universal scaling of the cumulants κ q (τ Q ) of the kink-pair number distribution P(n) as a function of the quench time. The experimental data (symbols) are compared with the scaling prediction (dashed lines) and the numerical data (solid lines) for the first three cumulants with q = 1, 2, 3. The universal scaling of the first cumulant κ 1 (τ Q ) = 〈n〉 is predicted by the KZM according to Eq. (4). Higher-order cumulants are also predicted to exhibit a universal scaling with the quench time. All cumulants of the experimental P(n) exhibit deviations from the universal scaling at long quench times consistent with a dephasing-induced anti-KZM behavior. Further deviations from the scaling behavior are observed at fast quench times. The range of quench times characterized by universal behavior is reduced for high-order cumulants as q increases. The error bars are much smaller than the size of the symbols used to depict the measured points. κ q (τ Q ) with q > 1 establishes the universal character of critical dynamics beyond the KZM. The later is explicitly verified for κ 2 (τ Q ), as shown in Fig. 4. Beyond the fast-quench deviations shared by the first cumulant, the power-law scaling of the variance of the number of kink pairs κ 2 (τ Q ) extends to all the larger values of the quench time explored, with barely noticeable deviations. By contrast, for the third cumulant κ 3 (τ Q ), the experimental data are already dominated by nonuniversal contributions both at short quench times, away from the scaling limit. For slow ramps, the third cumulant exhibits an onset of adiabatic dynamics due to finite-size effects around Jτ Q /ℏ = 10 2 . Finite-size effects are predicted to lead to a sharp suppression of the third cumulant. However, the experimental data shows that κ 3 starts to increase with τ Q . This is reminiscent of the anti-KZ behavior that has been reported in the literature for the first cumulant in the presence of heating sources 43,44 , and we attribute it to noise-induced dephasing of the trapped-ion qubit. The nature of these deviations is significantly more pronounced in high-order cumulants, reducing the regime of applicability of the universal scaling.

Discussion
In summary, using a trapped-ion quantum simulator we have measured the full distribution of topological defects in the quantum Ising chain, the paradigmatic model of quantum phase transitions, and characterized in detail its first three cumulants as a function of the quench rate. The statistics of the number of kink pairs have been shown to be described by the Poisson binomial distribution, with cumulants obeying a universal power law with the quench time, in which the phase transition is crossed. Our findings demonstrated that the scaling theory associated with a critical point rules the formation of topological defects beyond the scope of the KZM, which is restricted to the average number. Our work could be extended to probe systems with topological order, in which defect formation has been predicted to be anomalous 45 . We anticipate that the universal features of the full counting statistics of topological defects may be used in the error analysis of adiabatic quantum annealers, where the KZM already provides useful heuristics 10 .

Methods
A detailed description of analytical methods employed to obtain the results can be found in the Supplementary Information accompanying this work. In the following subsections, we present the experimental techniques employed. Experimental mapping at two-level system. Hamiltonian of a single qubit driven by a microwave field is given by Eq. (6). In the experiment, the Rabi frequency of the qubit simulator is set~20 kHz, which depends on the driven power of our microwave amplifier. The higher microwave power can shorten the Rabi time T R = 1/Ω R , and reduce the total operation time. To simulate the quenching process, one would like to vary g from −∞ to 0, an idealized evolution that needs infinite time, which cannot be realized in the experiment. However, to explore the universality associated with the crossing of the critical point one can initialize the system sufficiently deep in the paramagnetic phase. According to the KZM, it is sufficient to choose an initial value of the magnetic field such that the corresponding equilibrium relaxation time is much smaller than the time left until crossing the critical point. The system is then prepared out of the "frozen region" in the language of the adiabatic-impulse approximation 7 . For the quench time τ Q ≥ 1, the initial value of g = −5 is far out of the frozen time. We can simulate the TFQIM with an initial g = −5 and no excitations, by preparing the initial state of the qubit in the ground state ofĤ s k;i ¼Ĥ s k ðg ¼ À5Þ before the quenching process. The initial state can be derived by solving the eigenvectors of Eq. (6), which is where θ k;i ¼ Àarctan sin ka ð Þ 5 þ cos ka ð Þ . The scheme to detected the quantum critical dynamics of the one-dimensional TFQIM in each mode by using a single qubit is shown in Fig. 5. The whole process can be divided to three steps. Before the process, the qubit has been pumped to 0 j i state by using a 369 nm laser to excite transition S 1∕2 , F = 1 → P 1/2 , F = 1. In the first stage, the ion-trap qubit is prepared into the state ψ À k;i by a resonant microwave pulse. The second stage is the quench process; the Hamiltonian is time dependent and varies fromĤ s k;i Ĥ s k ðg ¼ À5Þ toĤ s k;f Ĥ s k ðg ¼ 0Þ. This is implemented by applying a chirped microwave pulse to the qubit. The chirped pulse is in the form of Þ is the pulse length. The third stage is to measure the excitation probability after the quench, which is to measure the occupation probability on the excited eigenstate of H s k;f . The excited eigenstate is ψ À k;f ¼ sin As the fluorescence detection on 171 Yb + can only discriminate 0 j i and 1 j i sates, we need to rotate the state ψ À k;f to 1 j i and then detect the bright state probability.
We thus use a qubit rotation and a fluorescence detection in this stage. The calculated waveform to set the AWG to perform the prerotation, quench and postrotation in the three stages will be discussed in the next section.
Driving waveform with an AWG. We consider the driving microwave from AWG is A cosðω c t þ ΦðtÞÞ, where A is the amplitude, ω c is the carrier frequency, and Φ(t) is an arbitrary phase function. The carrier frequency can be mixed up with some frequency ω 0 from a local oscillator by using a frequency mixer. The waveform at the output of the mixer is A cosðω c t þ ΦðtÞÞ cosðω o tÞ ¼ A=2ðcosððω o þ ω c Þt þ ΦðtÞÞ þ cosððω o À ω c Þt À ΦðtÞÞÞ: When the wave passes through the high-pass filter, the waveform is filtered as A=2 cosððω o þ ω c Þt þ ΦðtÞÞ. In the experiment, we choose the frequency ω o + ω c resonant with the qubit transition ω Q = ω o + ω c , so the magnetic field of the final driving waveform is BðtÞ ¼ B 0 cosðω Q t þ ΦðtÞÞ. The interaction between the magnetic field and spin is H I = −μ ⋅ B, where μ is the magnetic dipole of the spin qubit. The Hamiltonian can be further simplified after the rotating wave approximation in the interaction framê where the Rabi frequency Ω R ¼ À 0 h jμ Á B 1 j i=_. The phase function Φ(t) corresponds to the azimuthal angle in the Bloch sphere. Equation (6) can be also

expressed in interaction frame
so in the quench stage, the phase function is We rotate the state along a vector in the equatorial plane to prepare a state from 0 j i, or measure state to 1 j i, which means Φ(t) is constant in these two stages. The pulse length for the preparation and measurement stages is determined by the polar angle of the ground state at the beginning and end of the quench process, t = θ∕Ω R , respectively. The whole expression of Φ(t) in the three stages can be derived as ð0; t 1 Þ Ω R sin ka Ω R τ Q sin ka ðt À t 1 Þ 2 À 5 þ cos ka ð Þ ð t À t 1 Þ ; ðt 1 ; t 2 Þ ϕ f À π 2 ; ðt 2 ; t 3 Þ
Qubit preparation and detection error. There are some limitations to the preparation and measurement of the qubit. We measured the error through a preparation and detection experiment. First, we prepare the qubit in the 0 j i state by optical pumping method and detect the ion fluorescence. Ideally, no photon can be detected as the ion is in the dark state. However, there are dark counts of the photon detector as well as photons scattered from the environment. We repeat the process 10 5 times, and estimate the detected photon numbers. We also repeat the process 10 5 times by preparing the qubit in the bright state. Histograms for the photon number in the dark and bright states are shown in Fig. 6. In the experiment, the threshold is selected as 2: the state of the qubit is identified as the bright state when the detected photon number is ≥2. There is a bright error ϵ B for the dark state above the threshold, and a dark error ϵ D for the bright state below the threshold. The total error can be take as ϵ = (ϵ B + ϵ D )/2.

Data availability
The data that support the plots and other findings within this paper are available from the corresponding author upon reasonable request.
Received: 1 September 2019; Accepted: 29 January 2020; Fig. 6 Histogram for photon counts in state preparation and detection experiments. The distribution of photon counts is shown when the qubit state is prepared in 0 j i (dark state) and 1 j i (bright state).