Experimentally testing quantum critical dynamics beyond the Kibble-Zurek mechanism

TL;DR

This paper experimentally investigates the distribution of topological defects in a quantum phase transition, demonstrating universal scaling beyond the traditional Kibble-Zurek mechanism using a trapped ion quantum simulator.

Contribution

It extends the understanding of critical dynamics by experimentally verifying the full distribution of defects and their universal scaling laws beyond the mean, incorporating noise effects.

Findings

Kink pair distribution follows a Poisson binomial distribution.

Universal power-law scaling of cumulants with quench time.

Deviations due to noise-induced dephasing observed.

Abstract

We experimentally probe the distribution of kink pairs resulting from driving a one-dimensional quantum Ising chain through 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 the paradigmatic Kibble-Zurek mechanism, which accounts for the mean kink number, to characterize the full probability distribution of topological defects.

Equations

$$\hat{\mathcal{H}}=-J\sum _{m=1}^N\left({\hat{\sigma }}_m^z{\hat{\sigma }}_{m+1}^z+g{\hat{\sigma }}_m^x\right),$$

$$\hat{\mathcal{H}}=\sum _{k>0}{\hat{H}}_k=\sum _{k>0}{ϵ}_k\left(g\right)\left({\hat{\gamma }}_k^{†}{\hat{\gamma }}_k+{\hat{\gamma }}_{-k}^{†}{\hat{\gamma }}_{-k}-1\right),$$

$$g(t)=\frac{t}{{\tau }_Q}+g(0),$$

$$⟨n{⟩}_{\mathrm{KZM}}=\frac{N}{4\pi }\sqrt{\frac{\hslash }{2J{\tau }_Q}},$$

$$\tilde{P}(\theta )={\int }_{-\pi }^{\pi }d\theta P(n)e^{in\theta }=\prod _{k>0}\left[1+\left(e^{i\theta }-1\right)p_k\right],$$

$${\hat{H}}_k^{\mathrm{TLS}}=\frac{1}{2}\hslash \left({\Delta }_k(t){\hat{\sigma }}_z+{\Omega }_R{\hat{\sigma }}_x\right),$$

$$P(n)\simeq \frac{1}{\sqrt{6⟨n{⟩}_{\mathrm{KZM}}/\pi }}\exp \left[-\frac{{\pi }^2(n-⟨n{⟩}_{\mathrm{KZM}}{)}^2}{6⟨n{⟩}_{\mathrm{KZM}}}\right].$$

$$\hat{\mathcal{H}}=-J\sum _{m=1}^N({\hat{\sigma }}_m^z{\hat{\sigma }}_{m+1}^z+g{\hat{\sigma }}_m^x),$$

$${\hat{\sigma }}_m^z$$

$$\begin{array}{rl}\hat{\mathcal{H}}=& -J[\sum _{m=1}^{N-1}{\hat{c}}_m^{†}{\hat{c}}_{m+1}+{\hat{c}}_{m+1}^{†}{\hat{c}}_m+{\hat{c}}_m{\hat{c}}_{m+1}+{\hat{c}}_m^{†}{\hat{c}}_{m+1}^{†}+g\sum _{m=1}^N\left(1-2{\hat{c}}_m^{†}{\hat{c}}_m\right)\\ & +\hat{\Pi }\left({\hat{c}}_N^{†}{\hat{c}}_1+{\hat{c}}_N{\hat{c}}_1+{\hat{c}}_1^{†}{\hat{c}}_N+{\hat{c}}_1^{†}{\hat{c}}_N^{†}\right)].\end{array}$$

$${\hat{c}}_m=\frac{e^{-i\pi /4}}{\sqrt{N}}\sum _k\exp \left[im\left(ka\right)\right]{\hat{c}}_k,$$

$$k=\frac{\pi }{Na}\left(2j-1\right),j=-\frac{L}{2}+1,\dots ,\frac{L}{2}.$$

$${\hat{\mathcal{H}}}_{\mathrm{Even}}=-\sum _{0<k<\frac{\pi }{a}}\left(\begin{array}{cc}{\hat{c}}_k^{†}& {\hat{c}}_{-k}\end{array}\right){\hat{H}}_k\left(\begin{array}{c}{\hat{c}}_k\\ {\hat{c}}_{-k}^{†}\end{array}\right),$$

$${\hat{H}}_k\equiv 2J\left(g-\cos \left(ak\right)\right){\hat{\sigma }}^z+2J\sin \left(ak\right){\hat{\sigma }}^x.$$

$${\hat{\gamma }}_k=u_k{\hat{c}}_k+v_k{\hat{c}}_{-k}^{†},$$

$$\left\{{\hat{\gamma }}_k,{\hat{\gamma }}_{k^{\prime }}^{†}\right\}={\delta }_{k,k^{\prime }},\left\{{\hat{\gamma }}_k^{†},{\hat{\gamma }}_{k^{\prime }}^{†}\right\}=\left\{{\hat{\gamma }}_k,{\hat{\gamma }}_{k^{\prime }}\right\}=0.$$

$${\hat{\mathcal{H}}}_{\mathrm{Even}}=\sum _{0<k<\frac{\pi }{a}}{ϵ}_k\left(g\right)\left({\hat{\gamma }}_k^{†}{\hat{\gamma }}_k+{\hat{\gamma }}_{-k}^{†}{\hat{\gamma }}_{-k}-1\right),$$

$${ϵ}_0\equiv \frac{E_0}{NJ}=-\frac{1}{2\pi }{\int }_{-\pi }^{\pi }{ϵ}_k\left(g\right)dk.$$

$$\hat{\mathcal{N}}=\frac{1}{4}\sum _{m=1}^N(1-{\hat{\sigma }}_m^z{\hat{\sigma }}_{m+1}^z)=\sum _{k\ge 0}{\hat{\gamma }}_k^{†}{\hat{\gamma }}_k,$$

$$P(n)=\mathrm{tr}\left[\hat{\rho }\delta [\hat{\mathcal{N}}-n]\right],$$

$$\tilde{P}(\theta )$$

$$\log \tilde{P}(\theta )=\sum _{q=1}^{\infty }\frac{(i\theta {)}^q}{q!}{\kappa }_q.$$

$$\tilde{P}(\theta )=\prod _{k>0}\left[1+\left(e^{i\theta }-1\right)p_k\right],$$

$${\kappa }_1$$

$${\kappa }_2$$

$${\kappa }_3$$

$${\kappa }_1\equiv ⟨n⟩$$

$$=\frac{N}{2\pi }{\int }_0^{\pi }dk\exp \left(-\frac{1}{\hslash }2\pi J{\tau }_Q{k}^2\right).$$

$${\kappa }_1=\frac{N}{4\pi }\sqrt{\frac{\hslash }{2J{\tau }_Q}}\times \mathrm{erf}\left(\pi \sqrt{\frac{2\pi J{\tau }_Q}{\hslash }}\right),$$

$$⟨n{⟩}_{\mathrm{KZM}}=\frac{N}{4\pi }\sqrt{\frac{\hslash }{2J{\tau }_Q}},$$

$${\kappa }_2$$

$$=\frac{N}{2\pi }{\int }_0^{\pi }dk{e}^{-\frac{2\pi J{\tau }_Q}{\hslash }{k}^2}\left(1-e^{-\frac{2\pi J{\tau }_Q}{\hslash }{k}^2}\right).$$

Full text

*1**CAS Key Laboratory of Quantum Information, University of Science and Technology of China, Hefei 230026, China

2CAS Center For Excellence in Quantum Information and Quantum Physics, University of Science and Technology of China, Hefei 230026, China

3Donostia International Physics Center, E-20018 San Sebastián, Spain

4Departamento de Física, Universidad de los Andes, A.A. 4976, Bogotá D. C., Colombia

5IKERBASQUE, Basque Foundation for Science, E-48013 Bilbao, Spain

6Department of Physics, University of Massachusetts Boston, 100 Morrissey Boulevard, Boston, MA 02125

7Theoretical Division, Los Alamos National Laboratory, MS-B213, Los Alamos, NM 87545, USA*

Introduction

The 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 describes the nonadiabatic dynamics across a phase transition and predicts the formation of topological defects, such as vortices in superfluids and domain walls in spin systems [1]. Pioneering insights on the KZM were conceived in a cosmological setting [2] and applied to describe thermal continuous phase transitions [3, 4, 5]. The resulting KZM was latter extended to quantum phase transitions [6, 7, 8, 9]. Its central prediction is that the average total number of topological defects $\langle n^{T}\rangle$, formed when a system is driven through a critical point in a time scale $\tau_{Q}$, is given by a universal power-law $\langle n^{T}\rangle\sim\tau_{Q}^{-\beta}$. The power-law exponent $\beta=D\nu/(1+z\nu)$ is fixed by the dimensionality of the system $D$ and a combination of the equilibrium correlation-length and dynamic critical exponents, denoted by $\nu$ 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 one-dimensional 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.

The TFQIM is described by the Hamiltonian

[TABLE]

that we shall consider with periodic boundary conditions $\hat{\sigma}_{1}^{z}=\hat{\sigma}_{N+1}^{z}$ and $N$ even. Additionally, $\hat{\sigma}_{m}^{z}$ and $\hat{\sigma}_{m}^{x}$ 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|\gg 1$) and a doubly degenerated phase with ferromagnetic order ($|g|\ll 1$). There are two critical points at $g_{c}=\pm 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 seminal 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]

[TABLE]

where $\hat{\gamma}_{k}$ are quasiparticle operators, with $k$ labeling each mode and taking values $k=\frac{\pi}{N}\left(2m-1\right)$ with $m=-\frac{N}{2}+1,\ldots,\frac{N}{2}$. The energy $\epsilon_{k}$ of the $k$-th mode is $\epsilon_{k}\left(g\right)=2J\sqrt{\left(g-\cos k\right)^{2}+\sin^{2}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 operator $\hat{\mathcal{N}}~{}\equiv~{}\sum_{m=1}^{N}\left(\hat{1}-\hat{\sigma}_{m}^{z}\hat{\sigma}_{m+1}^{z}\right)/4$. The latter can be equivalently written as $\hat{\mathcal{N}}=\sum_{k>0}\hat{\gamma}_{k}^{\dagger}\hat{\gamma}_{k}$, where $\hat{\gamma}_{k}^{\dagger}\hat{\gamma}_{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 $\hat{\gamma}_{k}^{\dagger}\hat{\gamma}_{k}$ can be measured. To this end, we consider the quantum critical dynamics induced by a ramp of the magnetic field

[TABLE]

in a time scale $\tau_{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 [8, 32, 33]: the average number of kink pairs $\langle\hat{\mathcal{N}}\rangle=\langle n\rangle$ after the quench scales as

[TABLE]

in agreement with the universal power law $\langle n\rangle_{{\rm KZM}}\propto\tau_{Q}^{-\frac{\nu}{1+z\nu}}$ with critical exponents $\nu=z=1$ and one spatial dimension ($D=1$).

Results

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 observable $\hat{\mathcal{N}}$. 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 (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}=\langle\hat{\gamma}_{k}^{\dagger}\hat{\gamma}_{k}\rangle$. 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]

[TABLE]

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, $\langle n\rangle=\sum_{k>0}p_{k}$ and ${\rm Var}(n)=\sum_{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 $p_{k}=\exp\left(-\frac{\pi}{\hbar}J\tau_{Q}k^{2}\right)$ [8] (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}\mathrm{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 $\left|0\right\rangle\equiv\left|F=0,\,m_{F}=0\right\rangle$ and $\left|1\right\rangle\equiv\left|F=1,\,m_{F}=0\right\rangle$, as shown in Fig. 1b.

State preparation and manipulation: At zero static magnetic field, the splitting between $\left|0\right\rangle$ and $\left|1\right\rangle$ is 12.642812 GHz. We applied a static magnetic field of 4.66 G to define the quantization axis, which changes the $\left|0\right\rangle$ to $\left|1\right\rangle$ resonance frequency to 12.642819 GHz, and creates a 6.5 MHz Zeeman splittings for $\mathrm{{}^{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 around 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 around 12.242 GHz and 12.642 GHz, so a high pass filter is used to remove the 12.224 GHz wave. Finally, the wave around 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 $\left|0\right\rangle$ state, by applying a resonant light at 369nm 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 $\left|1\right\rangle$ is detected by fluorescence detection with another resonant light at 369nm, exciting ${S_{1/2}},\,F=1$ to ${P_{1/2}},\,F=0$. Fluorescence of the ion is collected by an objective with 0.4 numerical aperture (NA). A 935nm laser is used to prevent the state of the ion to jump to metastable states [41]. The initialization process can prepare the $\left|0\right\rangle$ state with fidelity $>99.9\%$. The total error associated with the state preparation and measurement is measured as 0.5%  [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 [8]. 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,

[TABLE]

where $\Omega_{R}=4J/\hbar$ and $\Delta_{k}(t)=4J[g(t)+\cos{k}]/(\hbar\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\tau_{Q}/\hbar$, 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/\Omega_{R}$ denoting the Rabi time. A typical operation process driven by the programmed microwave is illustrated on the Bloch sphere, see Supplementary Notes 3 and 4 for details. First, the qubit is prepared in the ground state of $\hat{H}_{k}^{\rm TLS}$ at $g(0)=-5$. Then, $g(t)$ is ramped from $-5$ to 0 with quenching rate $1/\tau_{Q}$. Finally, the ground state of the final $\hat{H}_{k}^{\rm TLS}$ is rotated to qubit $\left|0\right>$, so that the excited state is mapped to $\left|1\right>$, 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 (S11) 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]

[TABLE]

A comparison between $P(n)$ and Equation (7) is shown in Fig. 1. 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 non-universal 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,\dots$

To explore universal features in $P(n)$, we shall be concerned with the scaling of its cumulants $\kappa_{q}$ ($q=1,2,3\dots$) as a function of the quench time. We focus on the first three, theoretically derived in the scaling limit in [43]. The first one is set by the KZM estimate for the mean number of kink pairs $\kappa_{1}(\tau_{Q})=\langle n\rangle_{{\rm KZM}}$. The second one equals the variance and as shown in the SM is set by $\kappa_{2}(\tau_{Q})={\rm Var}(n)=(1-1/\sqrt{2})\langle n\rangle_{{\rm KZM}}$. Finally, the third cumulant is given by the third centered moment $\kappa_{3}(\tau_{Q})=\langle(n-\langle n\rangle)^{3}\rangle=(1-3/\sqrt{2}+2/\sqrt{3})\langle n\rangle_{{\rm KZM}}$. Indeed, all cumulants are predicted to be nonzero and proportional to $\langle n\rangle_{{\rm KZM}}$ [36], with the scaling of the first cumulant $\langle n\rangle_{{\rm KZM}}\propto\tau_{Q}^{-\frac{1}{2}}$ 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 Information, deviations from the power law occur at very fast quench times, satisfying $\tau_{Q}<\hbar/(\pi^{3}J)$. However, this regime is not probed in the experiment, as all data points are taken for larger values of $\tau_{Q}$. Within the parameters explored, deviations at fast quenches with $J\tau_{Q}/\hbar\sim 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 $\kappa_{q}(\tau_{Q})$ with $q>1$ establishes the universal character of critical dynamics beyond the KZM. The latter is explicitly verified for $\kappa_{2}(\tau_{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 $\kappa_{2}(\tau_{Q})$ extends to all the larger values of the quench time explored, with barely noticeable deviations. By contrast, for the third cumulant $\kappa_{3}(\tau_{Q})$ the experimental data is already dominated by non-universal 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\tau_{Q}/\hbar=10^{2}$. Finite-size effects are predicted to lead to a sharp suppression of the third cumulant. However, the experimental data shows that $\kappa_{3}$ starts to increase with $\tau_{Q}$. This is reminiscent of the anti-Kibble-Zurek behavior that has been reported in the literature for the first cumulant in the presence of heating sources [44, 45] 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 probed the full counting statistics of topological defects in the quantum Ising chain, the paradigmatic model of quantum phase transitions. The statistics of the number of kink pairs has 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 Kibble-Zurek mechanism, 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 [46]. 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 Kibble-Zurek mechanism already provides useful heuristics [10].

Methods

A detailed description of analytical methods and experimental techniques employed to obtain the results can be found in the Supplementary Information accompanying this work.

Data availability

The data that support the plots and other findings within this paper are available from the corresponding author upon reasonable request.

Acknowledgements

We acknowledge funding support from the National Key Research and Development Program of China (Nos. 2017YFA0304100, 2016YFA0302700), the National Natural Science Foundation of China (Nos. 61327901, 11774335, 11474268, 11734015, 11821404), Key Research Program of Frontier Sciences, CAS (No. QYZDY-SSW-SLH003), the Fundamental Research Funds for the Central Universities (Nos. WK2470000026, WK2470000018), Anhui Initiative in Quantum Information Technologies (Nos. AHY020100, AHY070000), and John Templeton Foundation.

Author contributions

A.d.C. proposed the project and developed the theoretical framework. J.-M.C. and Y.-F.H. took the experimental measurements. F. J. G.-R. developed numerical simulations and prepared the figures. All authors contributed to the analysis of the results and the writing of the manuscript. C.-F.L., G.-C.G. and A.d.C. supervised the project.

Additional information

Supplementary information accompanies this paper at url will be inserted by publisher.

Competing Interests

The authors declare that they have no competing financial interests.

Correspondence

Correspondence and requests for materials should be addressed to Y.-F. H. (email: hyf@ustc.edu.cn), C.-F. L. (email: cfli@ustc.edu.cn) and A. d. C. (email: adolfo.delcampo@dipc.org).

**—Supplementary Notes—

Experimentally testing quantum critical dynamics beyond the Kibble-Zurek mechanism

** Jin-Ming Cui1,2, Fernando Javier Gómez-Ruiz3,4,6, Yun-Feng Huang1,2,∗,

Chuan-Feng Li1,2,†, Guang-Can Guo1,2 & Adolfo del Campo3,5,6,7,‡

*1**CAS Key Laboratory of Quantum Information, University of Science and Technology of China, Hefei 230026, China

2CAS Center For Excellence in Quantum Information and Quantum Physics, University of Science and Technology of China, Hefei 230026, China

3Donostia International Physics Center, E-20018 San Sebastián, Spain

4Departamento de Física, Universidad de los Andes, A.A. 4976, Bogotá D. C., Colombia

5IKERBASQUE, Basque Foundation for Science, E-48013 Bilbao, Spain

6Department of Physics, University of Massachusetts Boston, 100 Morrissey Boulevard, Boston, MA 02125

7Theoretical Division, Los Alamos National Laboratory, MS-B213, Los Alamos, NM 87545, USA*

This PDF file includes:

• •

Supplementary Note 1. Mapping between spin and free-Fermion representations in the Ising model.

• •

Supplementary Note 2. Cumulants of the kink-pair number distribution $P(n)$.

• •

Supplementary Note 3. Experimental Setup.

• •

Supplementary Figure 1. Universal scaling of the cumulants $\kappa_{q}(\tau_{Q})$ of the kink-pair number distribution $P(n)$ as a function of the quench time.

• •

Supplementary Figure 2. Scheme to measure the excitation probability.

• •

Supplementary Figure 3. Histogram for photon counts in state preparation and detection experiments.

Supplementary Note 1 Mapping between spin and free-Fermion representations in the Ising

model

We consider the one-dimensional quantum Ising model in a transverse magnetic field $g$. The TFQIM Hamiltonian is given by

[TABLE]

where $J$ is the hopping parameter and $\hat{\sigma}_{m}^{\alpha}$ are the Pauli matrix along the direction $\alpha=x,z$ at site $m$. Here we assume for convenience that $N$ is even and impose periodic boundary conditions $\hat{\sigma}_{N+1}^{z}=\hat{\sigma}_{1}^{z}$. The Hamiltonian is mapped to a free fermion model through the Jordan-Wigner transformation [1]

[TABLE]

Inserting  (S2) into the Ising Hamiltonian Eq. (S1) we find that the Hamiltonian is quadratic in the Fermi operators

[TABLE]

This Hamiltonian includes terms such as $\hat{c}_{m}\hat{c}_{m+1}$ or $\hat{c}_{m+1}^{\dagger}\hat{c}_{m}^{\dagger}$ that violate the fermion-number conservation. However, the parity operator $\hat{\Pi}=\left(-1\right)^{\sum_{m=1}^{N}\hat{c}_{m}^{\dagger}\hat{c}_{m}}$ is a constant of motion having the value $+1$ and $-1$. We expand $\hat{c}_{m}$ into the Fourier modes

[TABLE]

where the pseudomomenta $k$ take values

[TABLE]

Using the choice of even $N$ and considering anti-periodic boundary condition $\hat{c}_{N+1}=-\hat{c}_{N}$, the Hamiltonian maps to

[TABLE]

where $\hat{H}_{k}$ is defined by

[TABLE]

Next, we use the Bogoliubov transformation to rewrite the Hamiltonian in terms of a new set of fermionic operators $\hat{\gamma}_{k}$ whose number is conserved. These new operators are defined via a unitary transformation in terms of the pair $\hat{c}_{k}$, $\hat{c}_{-k}^{\dagger}$:

[TABLE]

where $u_{k}$, $v_{k}$ are numbers satisfying $u_{k}^{2}+v_{k}^{2}=1$, $u_{-k}=u_{k}$, and $v_{-k}=-v_{k}$. It can be checked that canonical fermion anti-commutation relation for the $\hat{c}_{k}$ imply that the same relations are also satisfied by $\hat{\gamma}_{k}$, that is,

[TABLE]

The Hamiltonian is diagonalized in terms of $\hat{\gamma}_{k}$,

[TABLE]

where $\epsilon_{k}\left(g\right)=2J\sqrt{\left(g-\cos\left(ak\right)\right)^{2}+\sin^{2}\left(ak\right)}$, with the ground state defined $\hat{\gamma}_{k}\left|\textbf{0}\right\rangle=0$. The ground state free energy reads

[TABLE]

Supplementary Note 2 Cumulants of the kink-pair number distribution $P(n)$

The kink-pair number distribution $P(n)$ describes the probability of observing $n$ pairs of kinks. The number of kink pairs is measured by the observable

[TABLE]

with eigenvalues $n=0,1,2,\dots$ The kink-pair number distribution is thus defined as the quantum expectation value of the projector operator $\delta[\hat{\mathcal{N}}-n]$ on the subspace with $n$ pairs of kinks

[TABLE]

where the expectation value is taken over the final state $\hat{\rho}$ upon completion of the quench. The Fourier transform of $P(n)$ is the characteristic function

[TABLE]

The cumulants $\kappa_{q}$ of the distribution $P(n)$ are defined via the expansion of the cumulant-generating function, which is the logarithm of the characteristic function,

[TABLE]

We next focus on the characterization of $P(n)$ for the state that results from the nonadiabatic crossing of the phase transition in the TFQIM and derive the exact expression for the first few cumulants $\{\kappa_{q}\}$. To this end, we note that the nonadiabatic crossing of the phase transition results in the production of quasi-particles (and kinks) in pairs. Specifically, due to the conservation of momentum, quasiparticles with wavectors $+k$ and $-k$ are excited jointly. The distribution of the number of pairs of kinks is a Poisson binomial distribution associated with $N/2$ Bernoulli trials each with probability $p_{k}$ for $k>0$. Its characteristic function reads

[TABLE]

with $k\in\{\pi/N,3\pi/N,\ldots,(N-1)\pi/N\}$ for an Ising chain with period boundary conditions and $p_{k}$ denoting the excitation probability in the $k$-mode. The Poisson binomial distribution is well studied and its cumulants are known to be of the form

[TABLE]

Their explicit evaluation can be performed using the Landau-Zener formula for $p_{k}$. In the continuum limit, this yields

[TABLE]

In terms of the error function, it is found that

[TABLE]

an expression that is valid even for moderately fast quenches, as long as described by the Landau-Zener formula but possibly away from the scaling limit. However, for large $\tau_{Q}>\hbar/(2\pi^{3}J)$, ${\rm erf}\left(\frac{\sqrt{\pi}}{2d}\right)$ approaches unity and the average number of kink pair simply reads

[TABLE]

in agreement with the KZM [2].

A similar derivation shows that the variance is given by

[TABLE]

This results in the exact expression

[TABLE]

that in the scaling limit reduces to

[TABLE]

To characterize the non-normal features of the distribution of kinks, we also provide the exact expression for the third cumulant $\kappa_{3}$, that equals the third centered moment $\mu_{3}=\langle(n-\langle n\rangle)^{3}\rangle$. In particular, we note that $\kappa_{3}$ is related to the skewness $\gamma_{1}$ of the distribution, as $\gamma_{1}=\kappa_{2}/\kappa_{2}^{3/2}$. In the continuum, it reads

[TABLE]

This expression is simplified in the scaling limit, performing a power series expansion with $\tau_{Q}\rightarrow\infty$. To leading order, we obtain the expression quoted in the main text

[TABLE]

Note that the scaling expressions Eqs. (S18), (S20) and (S22) are obtained from the more general ones, Eqs. (S17), (S19) and (S21), whenever the error functions can be replaced by unity. The error function increases monotonically and swiftly from zero value to unity as the argument is increased from zero value. Indeed, it saturates at unity for fairly moderate values of the argument, e.g., ${\rm erf}(2)=0.995$. By imposing $\sqrt{\frac{2q\pi^{3}J\tau_{Q}}{\hbar}}>2$, we predict the onset of the scaling limit for $\kappa_{q}$ at

[TABLE]

This condition is satisfied in all the experimental realizations we have performed. As a result, it is justified to consider the scaling limit, in which it is actually possible to derived a closed form expression of the characteristic function [3]

[TABLE]

in terms of the polylogarithmic function ${\rm Li}_{3/2}(x)=\sum_{p=1}^{\infty}x^{p}/p^{3/2}$.

From this expression, it follows that all cumulants are actually nonzero and proportional to $\langle n\rangle_{{\rm KZM}}$. The kink-pair number distribution is manifestly non-normal. The Gaussian (normal) approximation discussed in the main text follows from neglecting all $\kappa_{q}$ with $q>2$. One can also approximate $1-1/\sqrt{2}\approx 3/\pi^{2}$, so that

[TABLE]

and thus

[TABLE]

The experimental data shows deviations from the scaling regime at short and long quench times that become more and more pronounced with increasing order $q$ of the cumulant $\kappa_{q}$. Deviations at long times are attributed to dephasing-induced heating and result in an anti-KZM behavior whereby the value of $\kappa_{q}$ increases with $\tau_{Q}$. Deviations at short times are of a different nature: For $J\tau_{Q}\sim 1$ deviations from the Landau-Zener formula occur, as shown in Fig. 2 of the main text.

For the sake of clarity we elaborate on the connection between the distribution of kink pairs and the distribution of the total number of kinks. To account for the latter, we note that its characteristic function is simply obtained as

[TABLE]

where the factor of $2\theta$ (instead of the $\theta$) accounts for the creation of kinks in pairs. Cumulants $\kappa_{q}^{T}$ of the total kink number distribution are simply related to cumulants $\kappa_{q}$ of the distribution of pairs by the identity

[TABLE]

Supplementary Note 3 Experimental Setup

The experimental task is to simulate the Hamiltonian of the $k$th mode in Eq. (S7) using the trapped ion qubit. According the the time-dependent Schrödinger equation, if we scale the time of evolution as $t^{\prime}=t\sin\left(ka\right)$, the Hamiltonian will be change to

[TABLE]

Eq. (S29) can be taken as a Hamiltonian of single qubit driven by a microwave field

[TABLE]

where

[TABLE]

is the detuning between the microwave frequency $\omega_{B}$ and the qubit transition frequency $\omega_{Q}$, and where $\Omega_{R}=4J/\hbar$ is the Rabi frequency under resonant microwave. In the experiment, the Rabi frequency of the qubit simulator is set around 20 kHz, which depends on the driven power of our microwave amplifier. Higher microwave power can shorten the Rabi time $T_{R}=1/\Omega_{R}$, and reduce the total operation time. To simulate the quench process, one would like to vary $g$ from $-\infty$ to 0, an idealized evolution that needs infinite time, which can not be realized in 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 [8]. For the quench time $\tau_{Q}\geq 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 $\hat{H}_{k,i}^{s}=\hat{H}_{k}^{s}(g=-5)$ before the quench process. The initial state can be derived by solving the eigenvectors of Eq. (S30), which is

[TABLE]

where $\theta_{k,i}=-\arctan\frac{\sin\left(ka\right)}{5+\cos\left(ka\right)}$ .

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. 2. The whole process can be divided to three steps. Before the process, the qubit has been pumped to $\left|0\right\rangle$ state by using a 369nm laser to excite transition ${S_{1/2}},\,F=1\rightarrow{P_{1/2}},\,F=1$. In the first stage, the ion-trap qubit is prepared into the state $\left|\psi_{-}\right\rangle_{k,i}$ by a resonant microwave pulse. The second stage is the quench process; the Hamiltonian is time dependent and varies from $\hat{H}_{k,i}^{s}\equiv\hat{H}_{k}^{s}(g=-5)$ to $\hat{H}_{k,f}^{s}\equiv\hat{H}_{k}^{s}(g=0)$. This is implemented by driving a chirped pulse into the qubit. The chirped pulse is in the form of Eq. (S31) with $g(t)=5(t/T_{p}-1)$, where $T_{p}=5J\tau_{Q}\sin(ka)/\left(\hbar\Omega_{R}\right)$ 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 $\hat{H}_{k,f}^{s}$. The excited eigenstate is $\left|\psi_{-}\right\rangle_{k,f}=\sin\theta_{k,f}\left|0\right\rangle+\cos\theta_{k,f}\left|1\right\rangle$, where $\theta_{k,f}=-ka$. As the fluorescence detection on ${}^{171}\mathrm{Yb}^{+}$ can only discriminate $\left|0\right\rangle$ and $\left|1\right\rangle$ sates, we need to rotate the state $\left|\psi_{-}\right\rangle_{k,f}$ to $\left|1\right\rangle$ 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 pre-rotation, quench and post-rotation in the three stages will be discussed in the next section.

Driving Waveform with an Arbitrary Wave Generator:

We consider the driving microwave from AWG is $A\cos(\omega_{c}t+\varPhi(t))$, where $A$ is the amplitude, $\omega_{c}$ is the carrier frequency and $\varPhi(t)$ is an arbitrary phase function. The carrier frequency can be mixed up with some frequency $\omega_{0}$ from a local oscillator (LO) by using a frequency mixer. The waveform at the output of the mixer is

[TABLE]

When the wave passes through the high pass filter, the wave form is filtered as $A/2\cos((\omega_{o}+\omega_{c})t+\varPhi(t))$. In the experiment, we choose the frequency $\omega_{o}+\omega_{c}$ resonant with the qubit transition $\omega_{Q}=\omega_{o}+\omega_{c}$, so the magnetic field of the final driving waveform is $\mathbf{B}(t)=\mathbf{B_{0}}\cos(\omega_{Q}t+\varPhi(t))$. The interaction between the magnetic field and spin is $H_{I}=-\mathbf{\mu}\cdot\mathbf{B}$, where $\mu$ is the magnetic dipole of the spin qubit. The Hamiltonian can be further simplified after the rotating wave approximation in the interaction frame

[TABLE]

where the Rabi frequency $\Omega_{R}=-\left\langle 0\right|\mathbf{\mu}\cdot\mathbf{B}\left|1\right>/\hbar$. The phase function $\varPhi(t)$ corresponds to the azimuthal angle in the Bloch sphere. Equation (S30) can be also expressed in interaction frame

[TABLE]

so in the quench stage, the phase function is

[TABLE]

We rotate the state along a vector in the equatorial plane to prepare a state from $\left|0\right>$, or measure state to $\left|1\right>$, which means $\varPhi(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=\theta/\Omega_{R}$, respectively. The whole expression of $\varPhi(t)$ in the three stages can be derived as

[TABLE]

where $t_{1}=(2\pi-\theta_{k,i})/\Omega_{R}$, $t_{2}=t_{1}+T_{p}$, $t_{3}=t_{2}+(2\pi+\theta_{k,f})/\Omega_{R}$ and $\phi_{f}=-5\tau_{Q}(2.5+\cos ka)$.

Qubit preparation and detection error:

There are some limitations in the preparation and measurement of the qubit. We measured the error through a preparation and detection experiment. First, we prepare the qubit in the $\left|0\right>$ state by optical pumping method, and detect the ion fluorescence. Ideally, no photon can be detected as the ion is in the 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 1000 times by preparing the qubit in the bright state. Histograms for the photon number in the dark and bright states are shown in the Supplementary Figure. 3. 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 greater or equal to 2. Obviously, there is a bright error $\epsilon_{B}$ for the dark state above the threshold, and a dark error $\epsilon_{D}$ for the bright state below the threshold. The total error can be take as $\epsilon=(\epsilon_{B}+\epsilon_{D})/2$.

Supplementary Notes: References

• [1] E. Barouch and B. M. McCoy, Phys. Rev. A 3, 786 (1971).
• [2] J. Dziarmaga, Phys. Rev. Lett. 95, 245701 (2005).
• [3] A. del Campo, Phys. Rev. Lett. 121, 200601 (2018).