# Quantum Kibble-Zurek physics in long-range transverse-field Ising models

###### Abstract

We analyze the quantum phase transitions taking place in a one-dimensional transverse field Ising model with long-range couplings that decay algebraically with distance. We are interested in the Kibble-Zurek universal scaling laws emerging in non-equilibrium dynamics and in the potential for the unambiguous observation of such behavior in a realistic experimental setup based on trapped ions. To this end, we determine the phase diagram of the model and the critical exponents characterizing its quantum phase transitions by means of density-matrix renormalization group calculations and finite-size scaling theory, which allows us to obtain good estimates for different range of ferro- and antiferromagnetic interactions. Beyond critical equilibrium properties, we tackle a non-equilibrium scenario in which quantum Kibble-Zurek scaling laws may be retrieved. Here it is found that the predicted non-equilibrium universal behavior, i.e. the scaling laws as a function of the quench rate and critical exponents, can be observed in systems comprising an experimentally feasible number of spins. Finally, a scheme is introduced to simulate the algebraically decaying couplings accurately by means of a digital quantum simulation with trapped ions. Our results suggest that quantum Kibble-Zurek physics can be explored and observed in state-of-the-art experiments with trapped ions realizing long-range Ising models.

## I Introduction

The diversity of the macroscopic forms of matter has driven the interest in a comprehensive understanding of the underlying physical processes. At equilibrium, we distinguish different phases in which matter can organize itself, i.e., with different order, depending on the external conditions that predominate. Accordingly, each phase exhibits a specific qualitative behavior such as the solid, liquid and gas phase of matter, the paramagnetic and the ferromagnetic phases of a material or the degeneracy of ground states leading into spontaneous symmetry breaking, which can typically be quantified by an order parameter Huang (1987). Similarly fascinating are the transitions from one phase to another when the external conditions change, which happen at certain values of external parameters such as the temperature or the strength of an applied magnetic field. The scientific advances taken in the past decades have led to an in-depth classification of phase transitions. An important class are continuous phase transitions which exhibit a continuous first derivative of the free energy and a discontinuity in the second (or higher order) derivative. Notably, for quantum phase transitions (QPTs), which take place at zero temperature, this condition is applied to the ground state energy Sachdev (2011).

The Ising model plays an important role in understanding critical phenomena. Even though it was originally regarded as too simple to account for magnetism due to the absence of a classical phase transition in one dimension Ising (1925), it is nowadays considered a cornerstone in classical and quantum statistical mechanics Huang (1987); Sachdev (2011); Dutta et al. (2015). Indeed, the transverse field Ising model (TIM) is arguably the simplest spin model which exhibits a continuous QPT, and constitutes the paradigmatic example in this realm Sachdev (2011); Dutta et al. (2015).

An essential aspect of continuous phase transitions can be described within the powerful framework of universality: When approaching the critical point, the physical behavior of the system becomes universal, i.e., for example the typical length scale of correlations and the order parameter follow a power law determined by the critical exponents that characterize the phase transition. Such critical exponents, together with the dimensionality of the system define a universality class Stanley (1999). In this manner, different systems that exhibit a phase transition belonging to the same universality class feature the same physical behavior close to the critical point. Furthermore, at a continuous QPT the energy gap often closes in a universal way Sachdev (2011).

Interestingly, universal behavior at a continuous phase transition is also present in time-dependent processes far beyond the equilibrium condition. A framework which describes universality in a time-dependent scenario is the Kibble-Zurek mechanism (KZM) of defect formation Kibble (1976, 1980); Zurek (1985, 1996). A continuous phase transition typically occurs at a symmetry breaking phase transition. The subject of the KZM is the time evolution with the system initially in an equilibrium state of the symmetric phase, where, subsequently, by changing the external parameter in time across the critical point the system is forced to select a symmetry broken configuration. The KZM then predicts that, in the final state of the evolution, a variety of non-equilibrium quantities scale as a power of the rate at which the system traverses the critical point. Crucially, the exponents of the scaling relations are determined by the equilibrium critical exponents of the phase transition. Hence, the KZM reveals an intrinsic connection between the equilibrium properties and the dynamics of the system. While classical systems have originally been studied, it has later been shown that the KZM can also be applied to QPTs Zurek et al. (2005); Damski (2005); Dziarmaga (2005); Polkovnikov (2005); Kolodrubetz et al. (2012); Chandran et al. (2012); Polkovnikov and Gritsev (2008) (see del Campo and Zurek (2014) for a review).

In many cases, the scaling laws can be obtained by resorting to the adiabatic-impulse approximation. In this approximation, the state remains frozen when the external parameter is close to the critical point, i.e., it does not follow the evolution of the Hamiltonian. Outside of this stage of the evolution the system is changing adiabatically. Within this simplified picture, the transition between the two stages is supposed to be sharp at a certain value of the external parameter which sets the length scale of the correlations in the final state. The scaling laws predicted by the adiabatic-impulse approximation have been investigated theoretically in a number of classical settings Laguna and Zurek (1997, 1998); del Campo et al. (2010); Damski and Zurek (2010); Nigmatullin et al. (2016); Puebla et al. (2017a); Silvi et al. (2016) and found to be in good agreement with experimental results Pyka et al. (2013); Ulm et al. (2013); Monaco et al. (2002); Navon et al. (2015); Beugnon and Navon (2017). The quantum Kibble-Zurek (QKZ) scaling laws, on the other hand, have been studied in Zurek et al. (2005); Polkovnikov (2005); Dziarmaga (2005); Damski (2005); Silvi et al. (2016).

Similar to many other quantum many-body phenomena Hohenberg and Halperin (1977); Polkovnikov et al. (2011); Eisert et al. (2015); Heyl et al. (2013); Heyl (2018); Jurcevic et al. (2017); Zhang et al. (2017a); Bernien et al. (2017), the experimental progress of the last years opens up a new approach to study the QKZ mechanism (QKZM) in a controlled way by employing quantum simulators. However, the high degree of control and protection against noise that is required to experimentally realize the KZM in a quantum system represents a significant challenge. Experimental confirmations of the QKZM have only recently been achieved with Bose-Einstein condensates Clark et al. (2016); Anquez et al. (2016), using Rydberg atoms Keesling et al. (2019), or by directly simulating the Landau-Zener dynamics Xu et al. (2014); Cui et al. (2016); Gong et al. (2016). In addition, it is worth mentioning other aspects of universality arising in non-equilibrium dynamics, which have been recently identified in Refs. Erne et al. (2018); Prüfer et al. (2018); Eigen et al. (2018). Despite this progress, various aspects of non-equilibrium dynamics in isolated and open many-body quantum systems remain to be explored Polkovnikov et al. (2011); Eisert et al. (2015).

On the other hand, the transverse field Ising model with long range interactions has been realized with great success in various experiments with trapped ions, see, e.g., Kim et al. (2010); Smith et al. (2016); Islam et al. (2011); Friedenauer et al. (2008); Kim et al. (2009); Islam et al. (2013); Richerme et al. (2014); Cohen et al. (2015); Neyenhuis et al. (2017); Britton et al. (2012); Safavi-Naini et al. (2018), and has been the subject of recent theoretical studies of, both, ground-state properties Koffel et al. (2012); Fey and Schmidt (2016); Vodola et al. (2016); Sun (2017); Defenu et al. (2017); Zhu et al. (2018) and non-equilibrium aspects such as entanglement dynamics Hauke and Tagliacozzo (2013); Schachenmayer et al. (2013); Pappalardi et al. (2018), the quantum Kibble-Zurek scaling Jaschke et al. (2017) and other dynamical properties Gong and Duan (2013); Titum et al. ().

In this article, we analyze the scaling laws predicted by the QKZM in the one-dimensional long-range transverse field Ising model (LRTIM) with ferro- and antiferromagnetic couplings that decay algebraically with distance, and discuss the feasibility of observing QKZM in a state-of-art trapped ion setup. To this end, first, the phase diagram of the QPT and the equilibrium critical exponents are calculated using density-matrix renormalization group (DMRG) White (1992); Schollwöck (2005, 2011). Employing finite-size scaling theory Fisher and Barber (1972); Brankov (1996), we obtain good estimates of these quantities for different interaction ranges, which we corroborate by means of finite-size collapse. Further, we compare our results with previously published works, Refs. Koffel et al. (2012); Fey and Schmidt (2016); Vodola et al. (2016); Sun (2017); Defenu et al. (2017); Zhu et al. (2018); Jaschke et al. (2017). For the antiferromagnetic case, inconsistent results have been reported previously: Our findings support the hypothesis that the LRTIM remains within the universality class of the nearest-neighbor Ising model. On the other hand, for the ferromagnetic case, our results are in good agreement with previous studies. In the second part of the article, we numerically investigate the QKZM scaling laws for various quantities, such as the number of defects. Here, we employ the exact simulation of a system with up to 18 spins. For these system sizes, which are well within reach for current experiments with trapped ions, we find reasonably good agreement with the scaling predicted by the QKZM and the adiabatic-impulse approximation. Our results indicate therefore that QKZ physics can be readily accessed in state-of-the-art experiments. Finally, we discuss an efficient scheme to implement the time-evolution under the model with trapped ions using a digital quantum simulator that exploits the natural structure of the interactions that are available for this platform.

The article is organized as follows. In Sec. II we introduce the LRTIM and determine its phase diagram and equilibrium critical exponents. In Sec. III we first review the KZM of defect formation in the quantum regime, to then show that the predicted QKZ scaling laws can be observed in the LRTIM with experimentally feasible parameters. Finally, in Sec. IV we discuss the experimental realization of the LRTIM with trapped-ions using a digital quantum simulator approach, while the main conclusions are summarized in Sec. V.

## Ii Equilibrium properties of the LRTIM

We investigate the LRTIM of a one-dimensional system consisting of an even number of spin-1/2 particles. The Hamiltonian of the model can be written as (with )

(1) |

with the strength of the transverse field , and algebraically decaying interactions

(2) |

where () for ferromagnetic (antiferromagnetic) couplings and . Here, denote the Pauli matrices acting on site . In the following we will choose . By controlling the parameter , different long-range interaction mechanisms can be described by the model, such as dipole-dipole () or van der Waals () interactions. Further, the limiting cases yield two important (solvable) models: First, in the limit the nearest-neighbor TIM is retrieved – the paradigmatic example of a system undergoing a QPT Sachdev (2011); Dutta et al. (2015). Second, for we have that and therefore the system is fully connected known as the Lipkin-Meshkov-Glick model Lipkin et al. (1965). In this case, the dimensionality of the system is .

Independent of the exact form of the couplings, the model exhibits a symmetry associated with the invariance of the Hamiltonian Eq. (1) under spin-flip, i.e., since , where

(3) |

denotes the parity operator. Thus, the Hilbert space is split into two subspaces with even and odd parity, respectively. Throughout the following, we will work within the positive parity subspace as it contains the ground state of for sufficiently large .

The ground state properties of the model depend crucially on the sign of , which can be illustrated with the help of the nearest-neighbor TIM: For antiferromagnetic couplings and the ground state is a superposition of two staggered-magnetic ordered states and , where (recall that we consider even). In contrast, for ferromagnetic couplings the ground state is a superposition of the two fully-magnetic states and . For the ground state becomes paramagnetic, i.e., a fully-polarized state where we introduced .

We now set out to determine the phase diagram of the ferro- and the antiferromagnetic LRTIM and characterize the universality classes for different values of the parameter . To this end, we employ that in the thermodynamic limit and close to the critical magnetic field strength various quantities are expected to show universal behavior. That is, for sufficiently close to the correlation length diverges as Sachdev (2011)

(4) |

while the energy gap and the order parameter vanish as

(5) |

and (with )

(6) |

respectively, where are the critical exponents. Here the operator associated with the order parameter takes the form

(7) |

and

(8) |

for ferro- and antiferromagnetic interactions, respectively. Additionally, we determine the Schmidt gap De Chiara et al. (2012); Lepori et al. (2013) as an indicator of the QPT, which is expected to behave as

(9) |

close to . The Schmidt gap is defined as the difference between largest two Schmidt coefficients for a specific bi-partition of the system. Here, we will consider a bipartition of the chain at the center into two blocks of equal size . The Schmidt gap is related to the entanglement spectrum, that is, the eigenvalues of the reduced density matrix, which can exhibit critical behavior Li and Haldane (2008) (however, it can be misleading, see Chandran et al. (2014)). Note that the entanglement spectrum Dalmonte et al. (2018), Renyi entropies for pure states Daley et al. (2012); Islam et al. (2015) and mixed state entanglement measures Marty et al. (2016); Cramer et al. (2011); Marty et al. (2014) are experimentally accessible.

For the nearest-neighbor interactions the critical point and the critical exponents can be determined exactly: One finds that with ferro- or antiferromagnetic interactions (these quantities coincide for both cases), the critical value is given by , while the exponents are and De Chiara et al. (2012). In particular, the order parameter and the Schmidt gap exhibit identical critical behavior.

To obtain the critical point and the exponents for the LRTIM we calculate , and using DMRG. Recall that for any finite system with spins, the symmetry is not spontaneously broken and thus, constraining to the positive parity subspace leads into . We therefore calculate instead which reveals critical behavior, and for which is expected. Here, to simulate the long-range interactions, we approximate the algebraic decaying couplings by a sum of exponentials Murg et al. (2010), see Appendix A for further information. We consider system sizes of up to spins with maximum bond dimension , such that numerical convergence was attained (see Appendix A).

In order to extrapolate the results obtained with systems of finitely many particles to the thermodynamic limit, we rely on finite size scaling theory Fisher and Barber (1972); Brankov (1996). Accordingly, any quantity whose behavior in the thermodynamic limit and sufficiently close to is of the form

(10) |

(where depends on ) can – for any finite – be written as

(11) |

Here, is the finite-size scaling function associated with which depends solely on the scaling variable . The term represents sub-leading corrections to the scaling with .

We can use the relation Eq. (11) to determine the phase diagram and the critical exponents. To find the critical field strength , we rely on the Binder cumulant Binder (1981),

(12) |

with () for ferromagnetic (antiferromagnetic) couplings. For the Binder cumulant, the relation Eq. (11) implies that close to

(13) |

That is, becomes size-independent for . Sub-leading corrections can result in small deviations from Eq. (13). They are taken into account as follows: First, for two systems of sizes and we determine as a function of the product the value at which . Then, using the hypothesis that for parameters and , we fit and the critical point to the values of that we have obtained for several and . See Appendix B for further details.

We repeat these steps for different values of and ferro- and antiferromagnetic couplings to obtain the phase diagram: Fig. 1(a) shows as a function of which sets the boundary between the paramagnetic (PM) and the ferromagnetic (F) or antiferromagnetic (AF) phase, respectively.

We remark that for ferromagnetic couplings we consider as for smaller values higher order finite-size corrections become more prominent and hinder the computation of the critical point with our approach. This is reflected by a decreasing value of as . In contrast, for antiferromagnetic couplings, similar behavior is observed for , i.e., for much smaller values of , see also Appendix B.

The critical exponents and for ferro- and antiferromagnetic couplings are obtained from the scaling with of, respectively, the energy gap , the squared order parameter (with () for ferromagnetic (antiferromagnetic) couplings) and the Schmidt gap at the critical point using the relation Eq. (11). On the other hand, to calculate the critical exponent we rely on the following scaling Ferrenberg and Landau (1991); Pelissetto and Vicari (2002); Angelini et al. (2014)

(14) |

We refer to Appendix B for the derivation of the previous expression, Eq. (14), as well as for an example showing the scaling (see Fig. 7(d)). In Fig. 1(b) we show the results obtained numerically as a function of . We find that for ferromagnetic couplings all the critical exponents significantly change with . On the other hand, for antiferromagnetic couplings the critical exponents remain very close to the values of the nearest neighbor Ising universality class.

Further, we corroborate the validity of most of our results by means of finite-size collapse and verify the quality of each collapse employing a chi-squared test, see Appendix C for a detailed discussion.

Notably, previous works have reported on the phase diagram of the LRTIM Koffel et al. (2012); Sun (2017); Zhu et al. (2018); Jaschke et al. (2017), determined the product based on linked-cluster expansion Fey and Schmidt (2016) and calculated the critical exponents and via finite-size scaling Sun (2017); Zhu et al. (2018) and through renormalization group techniques Defenu et al. (2017) (see also Ref. Sperstad et al. (2012) for Monte-Carlo simulations of a spin chain interacting with a non-Ohmic bath and its connection with long-range models). The results on ferromagnetic interactions (including Refs. Fey and Schmidt (2016); Defenu et al. (2017); Zhu et al. (2018); Jaschke et al. (2017)) and the phase diagram of the antiferromagnetic interactions (determined in Refs. Koffel et al. (2012); Vodola et al. (2016); Sun (2017); Jaschke et al. (2017)) are in good agreement with our findings. On the other hand, the universality class of the QPT in the antiferromagnetic case is still under debate. In particular, the results based on DMRG reported in Ref. Koffel et al. (2012) suggests that the model belongs to the nearest neighbor universality class for while for the critical exponents vary continuously with . Similarly, in Ref. Fey and Schmidt (2016) it is obtained that for the quantity is close to 1 whereas it deviates from 1 for smaller (although not consistently with Ref. Koffel et al. (2012)). In contrast, the results presented in Ref. Sun (2017) strongly suggest that the nearest neighbor universality class holds for the antiferromagnetic model for any , which agrees with our results (see Fey et al. (2019) for a similar observation in two dimensions).

## Iii QKZM in the LRTIM

In order to study the QKZM Polkovnikov (2005); Zurek et al. (2005); Dziarmaga (2005); Damski (2005) for the LRTIM, we consider the time evolution of the system where the external parameter is changed linearly in time according to

(15) |

with the value of the external parameter of the initial configuration and the duration of the evolution. In particular, throughout the ensuing, we consider the transition from the paramagnetic phase to the ferro- or antiferromagnetic phase. The state of the system that we are interested in is then given by the solution of the Schrödinger equation

(16) |

with the initial state , where denotes the ground state of .

Since the system undergoes a continuous phase transition it features a diverging relaxation time at the critical point , that in the close vicinity of (and in the thermodynamic limit) reads

(17) |

where accounts for the microscopic details. As a consequence, the time evolution ceases to be adiabatic when it approaches the critical point and excitations in the instantaneous eigenbasis of the Hamiltonian are created. As predicted by the QKZM, the formation of these excitations follows scaling laws that are determined by the equilibrium critical exponents of the QPT and that are observable in the behavior of various quantities that depend on the time-evolved state.

As mentioned in the Introduction, in order to derive the scaling laws, we rely on the adiabatic-impulse approximation which divides the dynamics into two regimes. In the adiabatic regime the state is supposed to adapt instantaneously to the Hamiltonian. On the other hand, in the impulse regime, the state does not follow the changes of ; it remains frozen. In the course of the time-evolution, adiabaticity is lost when the relaxation time exceeds the timescale on which the external parameter is changing Zurek et al. (2005); del Campo and Zurek (2014). For the protocol Eq. (15), the timescale takes the form

(18) |

The transition between the two regimes can thus be estimated to occur at values of the external parameter that satisfy the condition . This yields the scaling

(19) |

for the location and width of the impulse regime. We emphasize that this transition is just a working hypothesis, allowing for a heuristic derivation of the scaling with . More rigorous scaling approaches can substantiate this derivation even in cases where the heuristic approach fails Nikoghosyan et al. (2016). The QKZM now predicts that the average number of defects created during the time evolution is determined by the correlation length at the boundary, that is,

(20) |

where and denote the length and dimensionality of the system, respectively, and where we used Eq. (19) and the relation Eq. (4). Here, the operator corresponds to the number excitations in the eigenbasis of the final Hamiltonian as noted above. This quantity may however be difficult to access experimentally for general long-range interactions. We therefore consider the average number of domains

(21) |

which may be practically accessible. Here, the plus and the minus sign refer to ferro- and antiferromagnetic interactions, respectively. Since the operator defined in Eq. (21) commutes with the Hamiltonian for any value of , we expect for the same scaling as for given in Eq. (20) Polkovnikov (2005); Deng et al. (2009). Note that for the nearest-neighbor TIM with open boundary conditions as considered here, the two quantities are related as . In general, the final state is in a superposition of different spin configurations where each configuration comprises a particular number of domains, and thus, can be observed with a certain probability. We will come back to this in Sec. III.2. In case of nearest-neighbor interactions, the domain sizes of the configurations that mainly contribute to the final state are limited by , in agreement with the scaling Eq. (20). However, this argument does in general not apply to long-range interactions. For example, for , the Hamiltonian is invariant under permutations of the spins and therefore the number of the domains is not directly related to the number of defects.

Further, we consider two additional quantities whose QKZ scaling can be predicted by similar arguments Polkovnikov (2005); De Grandi et al. (2010a, b); Grandi and Polkovnikov (2010): The excitation probability,

(22) |

and the residual energy,

(23) |

where . The expected QKZ scaling of these quantities upon crossing the QPT are given by

(24) | |||

(25) |

while right at the critical point, i.e. at such that , their scaling becomes

(26) | ||||

(27) |

We refer the reader to Refs. De Grandi et al. (2010a, b); Grandi and Polkovnikov (2010) for a derivation of such scaling laws based on the adiabatic perturbation theory and adiabatic-impulse approximation. Notably, the residual energy Eq. (23) belongs to a larger class of quantities of the form that exhibit a QKZ scaling law, where denotes an operator for which there is critical behavior (i.e., fulfills Eq. (10)). Here, the scaling is given by and , where is the equilibrium critical exponent of Deng et al. (2009); De Grandi et al. (2010b, a); Grandi and Polkovnikov (2010). It is worth stressing that for systems QKZ scaling has been found only for ramps ending at the critical point Hwang et al. (2015); Puebla et al. (2017b); Defenu et al. (2018).

Notably, in order to be able to observe the QKZ scaling laws, the duration (of the linear quench) is required to be in a certain range. If the duration of the quench is too short, the time-evolution is not adiabatic, even away from the critical point. On the other hand, for a finite system, remains in the ground state if is considerably smaller than the energy gap at the critical point. The time evolution is in this case is fully adiabatic to a good approximation, and thus the excitation probability and the residual energy are expected to scale as Polkovnikov (2005); Rigolin et al. (2008); De Grandi et al. (2010a). For more details on the influence that a finite system size has on the scaling, see Sec. III.3.

### iii.1 Verification of the adiabatic-impulse approximation

As discussed in the previous section, the Kibble-Zurek scalings are obtained from the adiabatic-impulse approximation which is a strong simplification. We therefore first investigate numerically the validity of this approximation for our model. To this end, in order to determine the transition point , we measure the loss of adiabaticity of the time-evolution by the instantaneous ground-state fidelity

(28) |

When the system enters in the impulse regime, is expected to decrease suddenly. We thus identify the instant at which drops below a certain threshold . We determine for different quench times to obtain .

For the LRTIM with and and, both, ferro- and antiferromagnetic couplings, the scaling of is shown in Fig. 3. For this example, we employ the protocol Eq. (15) with . In the figure, we further show an exponential fit and the theoretically predicted scaling with (with and from Sec. II). We find that the fitted exponent follows closely the theoretically predicted universal values: The ratio between the exponents employing the equilibrium and the dynamical approach is given by for ferromagnetic and for antiferromagnetic couplings.

For lower values of , however, finite-size effects have a more significant impact on the scaling laws. For example, for and antiferromagnetic couplings we find , i.e., the values are no longer compatible within our error estimation. For ferromagnetic couplings the finite-size effects are even more prominent, as commented in Sec. II. In these latter cases, the numerically determined value of strongly depends on the selected threshold .

### iii.2 QKZ scaling laws

Next, we analyze the scaling laws of the excitation probability and the residual energy , defined in Eqs. (22) and (23), respectively. At the time , these quantities are expected to follow the scaling given in Eqs. (26) and (27). Further we investigate the scaling of the number of domains after traversing the critical point, , where for . To compare the predicted exponents with the numerical data we again calculate the time evolution under the protocol Eq. (15). For different values of and , we consider quench times and determine the scaling law by fitting the numerical data to , with labeling the three quantities, namely, , and . In Fig. 4(a) we show the typical behavior of these quantities as a function of for antiferromagnetic couplings with . The exponential scaling is fitted in the range . The fitted exponents are found to be very close to the expected scaling, ruling out trivial adiabatic scaling. We can clearly distinguish the regime of sudden quenches () where the scaling deviates from the one predicted by the QKZM, and the quasi-adiabatic regime (). Not visible in the figure is the adiabatic regime which occurs for longer quench times, which is characterized by and the trivial scaling of the and as discussed in Sec. III.

In Fig. 4(b), top panel, we show the ratio between the fitted exponent and the QKZ scaling prediction, for , as a function of . The strongest deviation from the equilibrium exponent occurs for the residual energy for small . This suggests that finite-size effects are more pronounced in this quantity. In the lower panel of Fig. 4(b) we plot the actual values for and for ferro- and antiferromagnetic couplings. The solid dark lines and the gray shaded region correspond to and its error for the equilibrium critical exponents and obtained in Sec. II. Our results support previous studies of the QKZM for the antiferromagnetic case Jaschke et al. (2017). In contrast, for ferromagnetic couplings we find that the exponents depend strongly on , which is not in agreement with the results reported in Jaschke et al. (2017).

The presence of finite-size effects is furthermore visible in the form of the distribution of the number of defects or domains del Campo (2018). That is, let denote the probability of observing domains at . Then, it is expected that forms a normal distribution if the quench times admit the Kibble-Zurek scaling, while it is of the form of a Poisson binomial distribution for slower quenches where the scaling breaks down due to finite-size effects. This behavior can precisely be observed in our model exemplified for and in Fig. 4(c): While for and the normal distribution fits the numerical data very well, for the data displays an exponential decay .

In addition, it is worth mentioning that the fitted exponents do not change significantly if the initial state is given by the fully-polarized state instead of the ground state . This can be important to experimentally observe the scaling laws, since the fully-polarized state may be simpler to be prepared.

### iii.3 Non-equilibrium finite-size collapse

The quantum Kibble-Zurek predictions are based on singularities at the critical point that, strictly speaking, are only present in the thermodynamic limit. In any finite system, we rely on the fact that the quantities of interest approach the critical behavior in a sufficiently rapid and systematic manner. For equilibrium quantities a powerful tool to approach this problem is the finite-size scaling theory that we have employed in Sec. II. A similar strategy can be also applied to the non-equilibrium scenario Kolodrubetz et al. (2012); Acevedo et al. (2014); Francuz et al. (2016); Chandran et al. (2012); Puebla et al. (2017a, b).

According to the QKZM, the length scale relevant for the scaling of a quantity with is given by the correlation length at . For any finite system of length , however, the correlation length cannot exceed the size of the system. Since

(29) |

we thus introduce a non-equilibrium finite-size scaling function depending on two scaling variables Francuz et al. (2016); Chandran et al. (2012); Puebla et al. (2017a, b) (cf. Eq. (11))

(30) |

Here, we suppose that the quantity is intensive in order to obtain a size-independent scaling function with and . Note that, for fully adiabatic dynamics, , the equilibrium scaling function must be recovered, thus (see for example Francuz et al. (2016); Puebla et al. (2017b)). The scaling function is therefore a non-equilibrium generalization of , and it is expected to be valid when the loss of adiabaticity occurs close to the QPT Francuz et al. (2016); Puebla et al. (2017b); Acevedo et al. (2014). That is, Eq. (30) does not hold for sudden or too fast ramps . On the other hand, for a fixed value of , the size-independent function contains the QKZ scaling () and the trivial quadratic scaling (). Following a similar argument as in Eq. (29), resorting to the adiabatic-impulse approximation, the relaxation time close to the QPT scales as , while its maximum value for a finite system follows . Hence, when both are comparable, that is, when , one expects a crossover between QKZ and the trivial scaling regimes.

In order to illustrate this, we compute non-equilibrium quantities after a time under the protocol Eq. (15) such that and thus, . In this case, it follows . We show the finite-size collapse for and for two representative cases, namely for ferromagnetic and for the antiferromagnetic couplings, plotted in Fig. 5(a) and (b), respectively. It is worth noting that the good collapse of the data further supports that universality already plays a significant role in the dynamics of systems comprising a few number of spins.

## Iv Experimental realization using a digital quantum simulator

Ising couplings with tunable interactions can experimentally be implemented with trapped ions. The theoretical proposal Porras and Cirac (2004); Deng et al. (2005) was followed by experiments comprising two Friedenauer et al. (2008) and, more recently, several tens of individually controllable ions realizing complex many-body physics, see, e.g., Refs. Smith et al. (2016); Richerme et al. (2014); Neyenhuis et al. (2017); Lanyon et al. (2017); Friis et al. (2018). In these setups, the energy levels of each ion that encode the qubit states are connected via either an optical Jurcevic et al. (2014, 2017); Lanyon et al. (2017); Friis et al. (2018) or a microwave transition Smith et al. (2016); Zhang et al. (2017a, b); Neyenhuis et al. (2017); Olmschenk et al. (2007). A state dependent force is obtained by applying a bichromatic laser field with frequencies of opposite detunings from the qubit transition frequency. The couplings of the resulting phonon-mediated effective spin-spin interaction are of the form

(31) |

where is the mass of the ions, are the amplitudes of the -th normal mode and the associated normal-mode frequency James (1998). The Rabi frequency of the -th ion and the wave number are parameters of the laser field whose exact form depends on the setup that is used. Additionally, a transverse magnetic field can be generated via an asymmetric detuning of the laser frequencies (AC-Stark shift) or with an additional laser field in resonance with the two qubit states.

The range of the couplings can be adjusted by changing the laser detuning from infinite range when is tuned close to the center-of-mass mode frequency to dipole-dipole interactions when is sufficiently far from any phonon mode frequency. The detuning can also be set to a value between these two cases such that the couplings approximately follows an algebraically decaying function as in Eq. (2); but the deviation from an ideal algebraic decay is typically strong and the quality of the approximation is not controllable.

On the other hand, the couplings given in Eq. (2) are highly useful to study long-range interactions in a systematic way. It is thus important to simulate, e.g., the Ising model with these couplings with a high degree of precision, which can be achieved using a digital quantum simulation. Indeed, as we show in the following, the time-evolution under this model can be simulated with controllable precision by means of a sequence of quantum gates that are implementable with trapped ions and can be constructed by exploiting the form of the interactions . Remarkably, our scheme can be used to implement any type of Ising couplings. As an example, we focus here on the implementation of algebraically decaying couplings and leave the investigation of the general case and of the efficiency of the scheme compared to other implementations for future work. In particular, we consider the time-evolved state of the model with a time-dependent magnetic field strength as it can be used for adiabatic ground-state preparation or to verify the QKZM. To simulate the time-evolution with a digital quantum simulator we divide it into a sequence of unitary operations using the Magnus method Blanes et al. (2009),

(32) |

with , and denoting the initial state. Here, on every interval , the time-evolution is approximated by means of a time-independent Hamiltonian,

(33) |

This can further be expanded by the Lie-Trotter product formula, such that we can write

(34) |

for operators , , with which are further specified below. Notably, the errors introduced in the two approximations, Eqs. (32) and (34), can exactly be bounded and are in both cases of order .

At first sight, the expansion in Eq. (34) may only seem possible with a quantum circuit that employs two-body gates. However, we show now that it can be realized using a quantum circuit of depth at most that is a combination of single-site operations and a global time-evolution under the naturally occuring interactions Eq. (31). Our scheme is thus an example of a hybrid digital-analog simulation of a quantum system Lamata (2017); Lamata et al. (2018); Arrazola et al. (2016). To derive the desired gate operations, we define the matrix whose entries are given by

(35) |

Here, we can choose in order that is positive semi definite with smallest eigenvalue 0. The eigendecomposition of then takes the form with labeled in decreasing order. Note that the Ising model with couplings given by and are equivalent since the diagonal elements of have no influence on any physical properties of the model. Moreover, the Eq. (34) is satisfied for

(36) |

for . Then, as already mentioned, the time-evolution under can physically be implement with trapped ions using local spin flips and the evolution under the Ising model with couplings . To this end, we first apply a spin-flip to any site with by means of the unitary operator

(37) |

Since all entries in are strictly positive, we can assume that and hence . Then, we employ the evolution under the Ising model which depends on the detuning and the Rabi frequencies where the latters are different for the factors in the Trotter decomposition. We choose with for all normal mode frequencies , and

(38) |

for some intensity . Here, we remark that in order to avoid phonon excitations we require for all normal mode frequencies that , with the Lamb-Dicke parameter. The time-evolution with these parameters then takes the form with

(39) |

for . Here, denotes the couplings given in Eq. (31) with Rabi frequencies . Finally, to complete the -th trotterization step the operator is applied again. The overall protocol to approximate the time-evolution operator using trapped ions reads

(40) |

where for and . Note that the gates consist only of single-site operations that can be applied simultaneously, while addresses the ions globally. The quantum circuit to simulate the evolution under algebraically decaying couplings with is depicted in Fig. 6 for ions, although it can be applied to systems of any number of qubits.

## V Conclusions

In summary, we have analyzed the QPTs taking place in the one-dimensional Ising model with algebraically decaying long-range interactions. We have determined the phase diagram for ferromagnetic and antiferromagnetic couplings as a function of the parameter , which accounts for the algebraically decaying couplings. For different values of , we have computed the equilibrium critical exponents of the QPTs in order to characterize their universality class. These results are obtained performing detailed simulations, based on DMRG calculations and finite-size scaling theory. We found that, while the paramagnetic-ferromagnetic QPT changes universality class for , the QPT for the frustrated case remains in its short-range universality class for . Our results partially support the findings of previous works. In addition, we have determined the critical exponents for the order parameter and Schmidt gap, which obey the same trend as and . The correctness of the critical point and exponents are tested by means of finite-size collapse. Having obtained the phase diagram and critical exponents, we tackled the non-equilibrium dynamics resulting from traversing the QPT by linearly changing the magnetic field strength in time. We focused on the emerging scaling laws as a function of the quench rate, as predicted by the KZM of defect formation in a quantum system. To this end, we have first verified the adiabatic-impulse approximation and then analyzed the scaling in the average number of domains, the residual energy and the excitation probability for various . As we show, the fitted exponents are very close to the predicted ones. We again corroborate our findings by means of finite-size scaling for the non-equilibrium states. Finally, we have discussed the experimental realization with trapped ions and proposed a scalable scheme to efficiently simulate the dynamics using a digital quantum simulator.

Our findings show that the KZM of the LRTIM can be observed in system comprising only a few tens of spins. In view of the current technological advances, the paradigmatic QKZM represents a phenomenon that can be useful to benchmark quantum simulators as well as to unveil novel physical properties of many-body systems. In particular, our results suggest that QKZM can be explored in state-of-the-art trapped-ion settings which provide a fully controllable platform that naturally exhibits long-range interactions, whose properties may be hard to simulate on a classical device already for moderate system sizes.

###### Acknowledgements.

R. P. thanks G. De Chiara for enlightening discussions and acknowledges the support by the SFI-DfE Investigator Programme (grant 15/IA/2864). This work was supported by the ERC Synergy grant BioQ and the authors acknowledge support by the state of Baden-Württemberg through bwHPC and the German Research Foundation (DFG) through grant no INST 40/467-1 FUGG.## Appendix A Details of the DMRG calculations

In order to find the ground state we employ a standard single-site variational matrix-product state method Schollwöck (2005, 2011). For this, the Hamiltonian has to be expressed as a matrix-product operator (MPO) with a bond dimension that is sufficiently small to store it. One approach to obtain a suitable MPO representation is to expand the couplings in terms of a sum of exponentials Murg et al. (2010) which we review in the following. That is, formally, one seeks for parameters and with

(41) |

In order to find these parameters, we can first consider the economy-sized QR-decomposition , where is an matrix whose entries are given by . The parameters are choosen to be the eigenvalues of , where () is the matrix containing the first (last) rows of . Here, denotes the pseudo-inverse of a matrix . The coefficients can then be found by a standard least-square optimization. In the numerical simulations, we utilize terms which results in a approximation of the couplings of the order or better.

Also, it can readily be shown Schollwöck (2011) that an operator of the form

(42) |

can be expressed as an MPO with bond dimension 3. Further, the sum of operators Eq. (42) exhibits and MPO representation of bond dimension at most .

The convergence of the variational algorithm is verified by monitoring the truncation error and the smallest value of the entanglement spectrum upon a bi-partition of size , denoted . In particular, we choose the bond dimension such that for . For , our maximum bond dimension of provides .

## Appendix B Additional information to the derivation of the ground state properties

It follows a description of the method employed to locate the critical point and determine the critical exponents. As commented in the main text, we resort to the Binder cumulant with , Eq. (12), to locate the critical point for various values. For two different system sizes and , the Binder cumulants intersect at , thus providing an estimate of the actual . An example is shown in Fig. 7(a). As , the intersection will occur closer to , and hence, by fitting these estimates to Angelini et al. (2014)

(43) |

we obtain precise results for , see Fig. 7(b). It is worth stressing that the exponent gives account for higher-order corrections (sub-leading finite-size scaling). In Fig. 7(c) we show the fitted parameters, which indicate that the ferromagnetic LRTIM exhibits stronger sub-leading corrections than its antiferromagnetic counterpart for . The fits for however become unstable, indicating that larger system sizes are required to obtain reliable estimates. At the critical point , we can obtain the critical exponents through finite size scaling. From Eq. (11), up to sub-leading corrections, we expect the scaling at the QPT for a quantity as described in the main text. The critical exponent however must be determined in a different manner. For , we achieve this by computing the derivative of ,

(44) |

Here, have used the finite-size scaling function, i.e.,

(45) |

as we are interested in the case . At the critical point, it follows

(46) |

Hence, can be obtained from the following expression: