Critical Phenomena and Kibble-Zurek Scaling in the Long-Range Quantum Ising Chain

We investigate an extension of the quantum Ising model in one spatial dimension including long-range $1 / r^{\alpha}$ interactions in its statics and dynamics with possible applications from heteronuclear polar molecules in optical lattices to trapped ions described by two-state spin systems. We introduce the statics of the system via both numerical techniques with finite size and infinite size matrix product states and a theoretical approaches using a truncated Jordan-Wigner transformation for the ferromagnetic and antiferromagnetic case and show that finite size effects have a crucial role shifting the quantum critical point of the external field by fifteen percent between thirty-two and around five-hundred spins. We numerically study the Kibble-Zurek hypothesis in the long-range quantum Ising model with Matrix Product States. A linear quench of the external field through the quantum critical point yields a power-law scaling of the defect density as a function of the total quench time. For example, the increase of the defect density is slower for longer-range models and the critical exponent changes by twenty-five per cent. Our study emphasizes the importance of such long-range interactions in statics and dynamics that could point to similar phenomena in a different setup of dynamical systems or for other models.


Introduction
Common quantum many-body models, such as the quantum Ising model or Hubbard models, include only a nearest-neighbor interaction, while e.g. the Coulomb forces or the dipoledipole interactions have power law tails of the form 1/r α . This leads to the core question: "to what extent does the physics of these models changes when including decaying long-range interactions?" Physical models with tunable long-range interactions have been successfully realized in recent experiments due to advances in atomic and molecular physics. This especially includes experiments in heteronuclear polar molecules in optical lattices [1,2] and trapped ions [3]. Further examples for the realization of a long-range Ising model are the cold ion experiment in [4] realizing a two dimensional lattice and ultracold molecules with longrange dipole interactions, which have been successfully realized in recent experiments [5,6]. Furthermore, realizations in Rydberg atoms through long-range van der Waals interactions [7,8] or solid state physics, e.g. based on LiHoF 4 [9], exist and provide a large number of spins, e.g. in comparison to ultracold ion experiments at the current stage. In addition, magnetic dipoles can lead to the same long-range interactions which have been realized in recent experiments in terms of Chromium [10], Erbium [11], and Dysprosium [12], and further setups with long-range interactions tuned by optical cavities have been made [13], too. The list of experiments can be further extended with nitrogen vacancies in diamonds with a ferromagnetic long-range Ising interaction [14]. Recent advances in quantum annealing have fostered the mapping of problems onto a generalized quantum Ising model [15], which is a more general Hamiltonian than the long-range interaction. This is as well the case for quantum spin glasses [16]. These systems have in common that they can be described by two-state spin models with 1/r α interactions coupled to external transverse fields, including the long-range quantum Ising (LRQI) model. While the experimental systems range from implementations in one to three dimensions, we study a one-dimensional system, because our numerical techniques are best suited to this case. The recent work of several groups discussed new emerging phases due to long-range interactions [17] and different regimes that hold or violate the Lieb-Robinson bounds, limiting the speed of transferring information depending on the initial state [18].
For numerical studies of many-body systems, tensor network methods offer a variety of algorithms which can be tailored to certain problems [19]. The infinite Matrix Product States (iMPS) [20] on infinite translationally invariant lattices and the Matrix Product States (MPS) for finite systems are well suited for the ground state calculations of quantum systems. Dynamics can be simulated with the Time-Dependent Variational Principle (TDVP) [21], where other methods as Krylov or local Runge-Kutta were explored for comparison [22,23]. The quantum Ising model is not only fundamental because of the various experimental implementations mentioned, but also since it has a similar classical model in the Ising model, originally introduced in 1925 [24] with analytic solutions in one and two dimensions. Moreover, the quantum model in d dimensions can be mapped onto the classical Ising model in d + 1 dimensions; the solution for the classical model in two dimensions is linked to our one dimensional quantum chain. Therefore, it is an ideal starting ground to explore long-range interactions before including them into other models and has always been the toy model to explore a huge variety of methods before using them on other models since it is the description for a vast majority of physical problems. This includes methods such as the exact spectrum, perturbation theory, or the mapping of the ferromagnetic case to the corresponding classical Ising model in higher dimensions without using numerical methods. It is applicable to many problems that can be mapped into two-level systems and therefore accessible for researchers from fields as different as solid state physics and atomic, molecular, optical physics just mentioned before.
While the experimental systems can be in higher dimension, we consider one dimensional spin chains with long-range interactions subjected to the external magnetic field in its transverse direction. From there, we lay out the discussions of static and dynamic simulations via tensor network methods for models with long-range interactions. Numerical algorithms based on tensor networks perform optimally in one dimension and have limitations in higher dimension, although they have been proposed for two dimensional systems [25]. However, the Kibble-Zurek mechanism yields a fascinating topic for the long-range quantum Ising model in real time evolution; for example, quenching across the critical point could be used for state preparation whenever it is easier to prepare the ground state in one phase.
As for the usual quantum Ising (QI) chain, that is integrable, both static properties such as the phase diagram, critical exponents or correlation lengths, and dynamic properties have been studied at length. The latter includes the Kibble-Zurek mechanism (KZM) being applied to show that the density of defects in the ferromagnetic chain follows a scaling law for linear quenches [26,27,28], i.e. the density of defects increases for faster quenches across the quantum critical point. Those properties of the Ising model have always paved the road to study other models. Then, there appears a natural question: how do these long-range interactions affect the well-known results in the quantum Ising model? These answers are relevant to ultracold molecules with characteristic long-range dipole-dipole interactions or Rydberg atoms with long-range interactions. Therefore, we present the reader new results on the Kibble-Zurek scaling with the comparison of the static results using three different approaches. Our analytical approach uses an approximation which truncates phase terms of the Jordan-Wigner transformation and the two numerical algorithms cover infinite and finite systems. With those comparisons we see finite size effects as the lower critical external fields for decreasing system sizes; dynamic simulations are based on these results. We introduce the dynamics with experimentally measurable observables such as the local magnetization in Fig. 1 and compare the scaling of defects with the Kibble-Zurek hypothesis. The Kibble-Zurek mechanism posits that for decreasing quench times, the defect density increases less for the ferromagnetic case in regard to the nearest-neighbor quantum Ising model. In contrast, the number of defects stays at the same rate for the antiferromagnetic model for faster quenches in comparison to the quantum Ising model. With the exponent of the Kibble-Zurek mechanism changing about 25% in the ferromagnetic case, this should be an observable effect. Apart from the prediction of defects when quenching across the quantum critical point for preparing states it might even allow for characterization of long-range effects according to the critical exponent.
Throughout the paper we concentrate on both the antiferromagnetic and ferromagnetic case. Moreover, we address the challenges with those studies, e.g. simulations breaking down in certain regimes of the long-range quantum Ising. We use the Open Source Matrix Product States package (OSMPS) [22] freely available from [29]. As a secondary goal, we provide with those practical calculations and discussions on the numerical results a basis for future researchers to interpret their results on statics and dynamics gained with OSMPS or any other tensor network method library. The models suitable for further studies range from more advanced spin models of XY Z-type to all kind of Hubbard models.
While static aspects such as the phase diagram of the long-range quantum Ising model have been studied in the past [30], the antiferromagnetic case has been investigated recently by Koffel [31] and extended in [32] to the correlations. We add to these works the aspects of infinite size MPS algorithms. A study in the thermodynamic limit with Linked Cluster Expansions for both the ferromagnetic and antiferromagnetic case has been presented in [33]. The dynamics of the long-range quantum Ising chain has lately attracted a huge interest. The entanglement growth of a bipartition during a quench was described in [34]. The speed of information spreading and the building of correlations was studied in [35,36,37,18,38] and the equilibration time of long-range quantum spin models in [39]. We not only extend the picture of dynamics in the long-range interaction in the Ising model to the Kibble-Zurek hypothesis, but provide a detailed path how such aspects can be studied numerically keeping a side-by-side view of the antiferromagnetic and ferromagnetic cases.
The structure of the paper is as follows: after introducing the Hamiltonian of the longrange quantum Ising model in Sec. 2, Sec. 3 presents the theoretical approximation for the phase boundary. The actual results are then discussed in Sec. 4 containing the phase boundaries for all different methods used, the observables during the dynamics, and the Kibble-Zurek scaling for the long-range quantum Ising model. We finish with the discussion of our results in the conclusion. The description of the numerical methods follows in the Appendix A and consists of the numerical algorithms used and the studies on their convergence. Further aspects on the convergence can be found in Appendix B.

Modeling Long-Range Interactions in the Quantum Ising Model
The system is modeled by the long-range quantum Ising Hamiltonian, where i, j are 1D lattice coordinates, J and the unitless h represent the strengths of the interaction and transverse field, respectively, and α takes some positive value parameterizing the long-range interactions. Examples for well-known models are Coulombic interactions with α = 1, dipole-dipole interactions decreasing with α = 3, and α = 6 for van der Waals models and the tail of the Lennard-Jones potential. In the following, we investigate phase diagrams both for the ferromagnetic (J > 0) and antiferromagnetic (J < 0) cases in Eq. (1). The long-range quantum Ising model in Eq. (1) reduces to the usual quantum Ising 1 σ x i , is a typical observable for spin models and realizable in cold ion experiments. The measurement during the quench from the paramagnetic phase to the to the ferromagnetic phase for different interaction strength governed by power law 1/|j − i| α displayed for α in steps of 0.2 is compared to the nearestneighbor quantum Ising (QI, dashed curve). The triangles mark spots with first maximum of the gradient disregarding numerical fluctuations. The magnetization in x-direction is vanishing equally for all α, since the quench is starting at h(t = 32, α) = 2h crit (α) and ending at h(t = τ, α) = 0. h crit is the value of the quantum critical point which differs for each α. We emphasize the unified dynamics of the quench independent of α for the symmetric setup around the critical point.
Hamiltonian in the limit of α → ∞, which exhibits a quantum phase transition at h = 1 both in ferro-and antiferromagnetic cases.
For h > 1 we have in both cases the paramagnetic phase where spins align with the external field in the limit h → ∞. The ferromagnetic phase appears for h < 1 where the ground state is described for h = 0 by a Z 2 symmetric ground state (|↑ . . . ↑ ± |↓ . . . ↓ ) / √ 2 subject to spontaneous symmetry breaking. The antiferromagnetic ground state has a staggered order in z-direction (|↑↓↑↓ . . . ± |↓↑↓↑ . . . ) / √ 2 referred to as Néel order. The Hamiltonian of the long-range quantum Ising in Eq. (1) builds the foundation of the following work and the standard quantum Ising from Eq. (2) is used frequently to check limits.
A common approach to quantum many-body systems are critical exponents. These universal scaling laws are independent of the system size and describe quantities as a function of the distance to the quantum critical point [40]. One typical example for such a quantity is the correlation length in the system. We relate to the critical exponents in the analysis of the Kibble-Zurek scaling, where [41] gives a review of the link between the critical exponents and the crossing of a quantum critical point in a quench. This setup is described by the critical exponent ν and the dynamical critical exponent z. The more detailed description is included in Appendix A.2.

Theoretical Boundary via a Truncated Jordan-Wigner Transformation
Having a general overview of the long-range quantum Ising model, we now derive a wellcontrolled analytic approximation to rely on. This is important to have a comparison to the numerical results introduced later. In general, the long-range quantum Ising problem is not exactly solvable like the nearest-neighbor case, but a truncation of the phase terms in the Jordan-Wigner transformation regains a long-range fermion model along the Kitaev model, which allows us to approximate the critical transverse field for quantum phase transitions, h crit . We show the derivation of the critical transverse fields with the use of the Jordan-Wigner transformation [42,40] and a truncation approximation. The derivation can be divided into two steps. First, we use the Jordan-Wigner transformation to map our spin Hamiltonian onto fermions. The Jordan-Wigner transformation is a non-local mapping from spins onto fermions, and fermions onto spins, where the non-locality of the transformation is necessary to obey the corresponding commutation relations. In the case of a nearest-neighbor spin Hamiltonian, the result of the Jordan-Wigner transformation is a nearest-neighbor fermion Hamiltonian where the non-local terms of the transformation cancel each other. In contrast, long-range spin Hamiltonians are mapped into long-range fermion Hamiltonians with k-body interactions of arbitrarily large k. After the truncation in the second step, we obtain a Kitaevtype Hamiltonian with long-range two-body interactions, which has been solved exactly in [43]. The Kitaev Hamiltonian is itself important as it describes spinless fermions with p-wave pairing interactions and displays topological states [44]. For a complete picture we show both steps to derive the theoretical boundary. In the following, we rewrite the long-range quantum Ising Hamiltonian in one dimension from Eq. (1) with a general site dependent coupling V (i, j): and keep a constant in front of the external field in order to have a unitless h. This general equation describes a variety of different models starting from a quantum spin glass with random V (i, j) originally proposed as classical model [45], over the completely connected case with V (i, j) = 1 for all i, j which is equivalent to a special case of the Lipkin-Meshkov-Glick model [46,47] ranging to chimera setup in the d-wave quantum annealing computers [48] with a limited number of interaction to other qubits. This Hamiltonian reduces to our long-range quantum Ising model with V (i, j) = −J/|i − j| α , while to the usual quantum Ising model with V (i, j) = −Jδ j−i,1 . In order to use the Jordan-Wigner transformation, we have to rotate our coordinate system around the y-axis mapping σ x → σ z and σ z → −σ x following the approach used in [40]. We obtain where the minus sign of the transformation cancels out appearing as a pair in the interaction term. The Jordan-Wigner transformation, which is exact for the usual quantum Ising model, is given by where c i , c † i are annihilation and creation operators for spinless fermions satisfying the anticommutation relations. Applying the Jordan-Wigner transformation, we can rewrite the Hamiltonian Eq. (4) with the fermionic operators, Note that for the nearest-neighbor potentials, like in the usual quantum Ising model, there are no operator products in terms of the index m, thus the Hamiltonian becomes quadratic in the fermionic operators, reducing to the well-known Jordan-Wigner transformed quantum Ising model. The approximation step replaces the operator products m (1−2c † i c i ) in the interaction term of the Hamiltonian running over the index m with an identity. This is equivalent with truncating fermionic operators up to quadratic order to obtain an analytically solvable form and the reason behind calling it a truncated Jordan-Wigner approach. We continue with This Hamiltonian corresponds to the long-range Kitaev model exactly solved in [43]. Our truncation approximation becomes exact only for the nearest-neighbor potentials. In general, it becomes better for ferromagnetic states, since they satisfy We perform a standard Fourier transform and Bogoliubov transform [40], which completes the diagonalization of Eq. (9) with the result, where the excitation spectrum q is given by with cosine and sine transforms of the potential where we replaced (i, j) with the distance r = j − i valid for the long-range quantum Ising model respectively. The excitation spectrum Eq. (11) could have a gapless Γ-point (q = 0) at the critical field and a gapless K-point (q = ±π) at though it depends on details of the potentials as to whether or not the spectrum exhibits the gapless points. For ferromagnetic long-range interaction, where ζ(α) is the Riemann Zeta function. In contrast, for antiferromagnetic long-range interaction, V (r) = |J|/r α , we have a gapless K-point at We remark that for the usual quantum Ising model with V (i, j) = −Jδ j−i,1 , we expect a gapless Γ-point at h Γ crit = 1 in the ferromagnetic case, and a gapless K-point at h Γ crit = 1 in the antiferromagnetic case.

Phase Diagram of the Long-Range Quantum Ising Chain
Our first result for the phase boundary is obtained analytically. In Sec. 3 we used the Jordan-Wigner transformation, which solves the nearest-neighbor quantum Ising chain exactly, for the long-range quantum Ising leading to an infinite series of fermionic operators. We truncate such series to obtain analytic estimates for the critical magnetic fields and recall the results from Eqs. (16) and (17) for the ferromagnetic and antiferromagnetic cases, respectively. We compare this to the numerical values of the critical field h crit discussed thoroughly in Appendix A.1. In Fig. 2 we show the phase boundaries for the different methods, where the results of the Kitaev model from Sec. 3 corresponds to the solid blue curves. We observe as a first trend that in the ferromagnetic long-range quantum Ising chain, as α decreases, that is, as the potential becomes more strongly long-range, the ferromagnetic phase becomes more favorable: the negative energy contribution from the ferromagnetic ordering increases in its absolute value. On the other hand, in the antiferromagnetic long-range quantum Ising chain, the antiferromagnetic phase becomes less favorable for smaller values of α because the long-range potential induces more frustration in the staggered ordering. Next, we perform numerical simulations to estimate the critical lines more accurately.
During the numerical simulations we follow two different paths. First we use the infinite Matrix Product State algorithm [20,22] to search for a translationally invariant ground state for a unit cell of q sites. We search for this kind of ground state for q = 2, although the Hamiltonian in Eq. (1) fulfills translational invariance in the thermodynamic limit for any q. The energy is minimized via variational methods and the strong finite size effects at the boundaries due to the long-range interactions are avoided. Second, we consider finite size systems with a number of sites L ∈ {32, 64, 128, 256, 512} and approximate the critical point in the thermodynamic limit L → ∞ via finite size scaling. On the one hand this is an appropriate check on the results. On the other hand, we use finite size systems for the dynamics later on and provide a prediction for finite size effects commonly present in quantum simulators. To identify the critical point we keep α and J fixed and iterate over the transverse field h. The von Neumann entropy as a function of the external field has a maximum located at the critical value of the external field h crit . For the nearest-neighbor model the critical point can be found based that the area law for entanglement is violated for gapless states [49], where the entanglement is measured with the von Neumann entropy. Although long-range models may in general violate this area law away from the critical point, we take the same approach knowing that the ground state obeys the area law for the two limiting cases h → 0 and h → ∞. The Ising model has a closing gap around the critical point and therefore the von Neumann entropy S can be used as indicator for the quantum critical point where λ i are the singular values squared of the Schmidt decomposition of the system and χ is the maximum bond dimension controlling the truncation of the Hilbert space. In the finite size systems, the eigenvalues of the reduced density matrices correspond to λ i when splitting into a bipartition. The entropy and the entanglement have a unique extremum at the phase transition which must coincident with the closing gap and can therefore be used to identify the critical point. An alternative approach would be to use the correlation length of the system or the maximal gradient in the magnetization, both with respect to the external field. We explore the critical behavior of the infinite chain (triangles) and the finite size scaling results (diamonds) in Fig. 2.
In conclusion, the relative difference between the iMPS and the theoretical estimate is bound by 0.37 (0.17) for the ferromagnetic (antiferromagnetic) model considering α ≥ 2. The relative difference for two arbitrary values x and x is defined as We compare the theoretical approximation (Theo.) based on a truncated Jordan Wigner transformation, the infinite MPS results (iMPS), and the finite size scaling (FSS) of the MPS results. The boundary between the gray and white shading is the critical point in the limit α → ∞ corresponding to the nearest-neighbor quantum Ising model. We observe that both cases react differently when increasing the long-range interactions. In the ferromagnetic case the paramagnetic phase decreases at cost of the ferromagnet, while the antiferromagnetic phase becomes smaller in the other case due to the concurrency in the staggered order. All three different methods show these trends, although they show in the regime α significant relative differences above 0.1.
The actual trend with the difference between iMPS and the theory curve vanishing in the limit α → ∞ reflects the analysis in Sec. 3. Further we state that the relative difference between iMPS and finite size scaling results are below 0.11 (0.04). The details on the finite size scaling are described in Appendix A.1.

Scaling of Defect Density for the Kibble-Zurek Hypothesis under a Quantum Quench
We now turn our focus to the dynamics of the long-range quantum Ising model. We analyze the quench through the quantum critical point by means of finite systems. We first analyze measurable observables during the quench and then determine the Kibble-Zurek scaling according to the final state, that is a power-law describing the generation of defects as a function of the quench time and a critical exponent. The linear quench for a time dependent external field is defined as and replaces the coupling h in Eq. (1). τ is the total time of the quench. We quench from When considering our observables during the quench, we stay close to experimentally realizable measurements. In Fig. 1 we presented the average magnetizations x This is e.g. measurable in cold ion experiments [4]. We choose to plot the magnetization in the x-direction which is measurable in the ferromagnetic and antiferromagnetic case, and can not average out in the simulations due to superpositions, e.g. the ferromagnetic ground state α → ∞, and h = 0: (| ↑ . . . ↑ ± | ↓ . . . ↓ )/ √ 2 has no net magnetization for a local measurement of σ z . We note in Fig. 1 for α ≈ 4 the steep gradient around t/τ ≈ 0.5, which corresponds to the critical value for the transverse field in the ferromagnetic case treated here. Therefore, we extract the first maximum of the gradient disregarding all minor fluctuations and plot them over α in Fig. 3. We see a similarity to the critical value of the transverse field from the static calculations for big τ , while for faster quenches the deviation from the static result increases. We present the analogous analysis for the antiferromagnetic case as well in Fig. 3. For faster quenches the maximal gradient is again below the static result leading to the assumption that this effect is due to a delay during which the system takes time to react. The antiferromagnetic case shows an additional effect as a function of α: longer-range interactions exhibit an extra delay in their maximum gradient. If we consider the antiferromagnetic quench for τ = 128 in Fig. 3 (a), the maximum of the gradient moves from h = h crit for α = 4 to values of h ≤ 0.75h crit for α < 2. One suggestion to explain this is the concurrent spin configuration in the staggered order in the antiferromagnet which grows for smaller α and might cause this extra delay to establish an antiferromagnetic state. The results on this analysis vary if the quench scheme is not symmetric around the critical point h crit .
Next, we present results for a non-global spin measurement. Such a measurement of a single spin in an optical lattice is technically possible [50] and can be realized in cold ion experiments, too. In Fig. 4 we show the nearest-neighbor spin-spin correlation in the z-direction for the antiferromagnetic case as a function of time and site in the system. As mentioned above, local measurements of σ z average out due to symmetry apart from numerical fluctuations. Correlations grow during the quench, and finite system boundary effects are visible at both ends of the chain. The total quench time is τ = 32, the system size L = 128, and α = 3 in this example, corresponding to dipole-dipole interactions. Now we turn to the Kibble-Zurek hypothesis in the quenches already discussed above. The Kibble-Zurek mechanism [26,51,41] predicts a scaling of the defects in the quantum system as a function of the quench time τ . We recall that the scaling for the nearest-neighbor Ising model is [52]: where we use ρ for the defect density, and ν and z for the critical exponents. We have already considered a one-dimensional system in Eq. (22). Furthermore, Eq. (22) and the following analysis applies only to linear quenches defined over ∂h/∂t = const., but extensions for the Kibble-Zurek scaling to non-linear quenches are possible [53]. Plugging in the critical exponents for the nearest-neighbor quantum Ising model with ν = 1 and z = 1 we obtain ρ ∝ τ −0.5 . The exponent 1/2 is a checkpoint for the calculations in the nearest-neighbor limit α → ∞. The match between theory and experiment does depend on the experimental realization, but yields all over a power law for the defects according to the review in [41].
We are interested in how much the exponent changes when tuning α for the long-range interactions and how it deviates from the 1/2 in the nearest-neighbor case. We estimate the exponent µ in the proceeding analysis which is defined as Although the scaling can be used in more general cases, we restrict ourselves to the linear quenches crossing the quantum phase transition located at the critical value h crit analyzed in the previous Section 4. This quench scheme starting at twice the critical value was discussed in Eq.   [21] at χ = 95, 100. We expect that the critical exponent µ approaches 0.5 in the limits α → ∞ and L → ∞ resulting in the ferromagnetic nearest-neighbor quantum Ising model in the thermodynamic limit. We provide a brief overview of the finite size effects in nearestneighbor model in Appendix D to support this statement. The error bars are the standard deviation for µ of the fit for the different values of τ . Examples for the fit are shown in Fig. A4 in Appendix A.2. It does not contain the truncation or methodical error of the time evolution method. From the two curves with the different bond dimensions we deduce that certain simulations fail, resulting in a large standard deviation in the fitting procedure. We consider those points with large standard deviation as a failure of convergence since another set of simulations with a slightly different bond dimension fit into the trend of converged points and have an equally small standard deviation. The critical exponent grows for increasing values of α, meaning if we double the quench velocity for two systems with different long-range interactions, the number of defects increases more slowly in the system governed by longerrange interactions. But as we show in Appendix A.2, longer-range systems have initially a higher defect density as shown in Fig. A4.
However, these results just discussed for the ferromagnetic case are not true for the antiferromagnetic model. The critical exponent stays constant independent of the long-range interaction. The deviation from the nearest-neighbor case µ = 0.5 should be due to finite size effects. But we know from Fig. A4 that the density of defects is again higher for longer range interactions. The defect density of the antiferromagnetic model is defined analogously to the ferromagnetic model: a missing kink in the staggered order is one defect.

Conclusion and Open Questions
In this study we have established that the introduction of long-range interactions changes the physics of one-dimensional quantum spin systems non-trivially as observed in various aspects of the statics and the defect density and therefore merits further attention. This is not limited to the quantum Ising model, e.g. the Bose-Hubbard and many other models could be discussed under the influence of the long-range interactions such as dipole-dipole interactions. First we summarize the results from our study of the long-range quantum Ising model and then raise open problems associated with long-range quantum physics.
We have presented multiple aspects of long-range quantum Ising model, e.g. the Kibble-Zurek scaling, which lead to a better understanding of physical systems with long-range interactions, e.g. the scaling of generated defects during the quench through a quantum critical point. For the ferromagnetic case we find a slower increasing defect density in comparison to the nearest-neighbor limit when increasing the quench velocity. The defect density is higher in comparison to the nearest-neighbor model for the same quench time. Our simulations give a detailed picture of the critical value of the transverse field in the thermodynamic limit using infinite Matrix Product States along with a practical guidance how to resolve issues introduced by degeneracy in those systems. We also presented an alternative to the infinite size algorithm using finite size scaling methods to obtain the value of the critical external field in the thermodynamic limit. This uncovers the finite-size effects, which will certainly be present in experiments with small system sizes. The relative difference for the value of the critical point between L = 32 and L = 512 is approximately 15% for the antiferromagnetic case for α = 3. The same comparison for the ferromagnetic yields a relative difference of around 20%. The iMPS and finite size scaling phase boundaries are supported by an approximate analytical result using a truncated Jordan-Wigner transformation approach in combination with the Kitaev model. We introduced the dynamics of the long-range quantum Ising model through local and global observables during a quench from the paramagnetic phase into the ferromagnetic and antiferromagnetic phases crossing the critical point. We related those measurements to the actual quantum critical point through the maximum of the gradient of average global magnetization in the x-direction and analyzed how this is affected by the total time of the quench. Finally, we evaluated the scaling of defects in the Kibble-Zurek scenario leading to the result of a slowly increasing (constant) defect density at the end for increasing quench velocity in comparison to the nearest-neighbor limit of the ferromagnetic (antiferromagnetic) long-range quantum Ising model. This is reflected in the exponent of the Kibble-Zurek scaling, e.g. in the ferromagnetic case we observe exponents in the range from µ = 0.45 to µ = 0.6 for power law exponents α > 1.5, where µ = 0.5 is the limit of the nearest-neighbor quantum Ising model in the thermodynamic limit.
Having concluded our study, we point out open problems that can be addressed in future research. As indicated by the list of the recent studies concerning the dynamics of the longrange quantum Ising model [34,36,37,18], static physics has been well explored [30,31,32], but near-equilibrium dynamics such as the Kibble-Zurek hypothesis remain largely unknown, let alone far-from-equilibrium dynamics. One aspect which could be explored in statics and dynamics is the extension to two or higher dimensional systems which we have excluded due to the limitation of MPS algorithms to one spatial dimension for studies of this kind. Twodimensional setups allow more options for quantum computation to apply gates [48]. In the case of the antiferromagnetic case triangular lattices can lead to the similar concurrence of the staggered order than for long-range interactions in the one-dimensional case. While twodimensional systems might be accessible through simulations with tensor network methods such as Projected Entangled Pair States (PEPS) [25], this route seems to be more accessible for experiments for any of the models mentioned before. Discrete Truncated Wigner methods are another tool to approach two-dimensional systems with numerical simulations [54]. A further path starting from here is the investigation of dynamical phase diagrams. Where recent studies treat e.g. binary pulses in regards to interaction and transverse field of the quantum Ising model [55], further possibilities considering time-dependent periodic long-range interactions arise. One study obtaining a dynamical phase diagram for the long-range quantum Ising chain in the infinite limit using iMPS is [56]. So called completely interacting models at α = 0 form another topic for research.
We have shown that the long-range quantum Ising model is a good starting point for any kind of long-range study. Other models may show new behaviors when including longrange interactions. We think here of the generalization of the quantum Ising in terms of the XY -model with a realization in Rydberg systems [57] or the XY Z-model as realized in polyatomic molecules [58]. Furthermore, this includes Hubbard models such as the Bose-Hubbard model, Fermi-Hubbard model as they appear e.g. for molecules. For macroscopic quantum tunneling in the Bose-Hubbard model [59] results might change if the distance of the long-range interaction is larger than the width of the barrier. Other problems suitable for extension to long-range interactions include the quantum rotor model and the bilinear biquadratic spin-one Hamiltonian [60]. of the existing system and variationally optimized in order to converge to the ground state. The previous iteration is included as the environment, where environment is understood in the sense that it exchanges quantum numbers with the system, and not in an system-bath open quantum system context. In the first step of the algorithm, the unit cell is optimized without any environment; in the second step the solution of the first step is considered as the environment. From that point on the environment grows subsequently with each step. The algorithm terminates if the newly introduced sites fulfill the orthogonality fidelity condition: where n indicates the iteration step and R is the right part of the system excluding the site introduces in the most recent step. Therefore, the density matrices ρ n−1 and ρ R n represent the same part of the system and the latter one has one more update included. When the overlap is sufficiently big, the update did not change the result any more. We refer to [20,22] for a detailed description. The phase boundary for the infinite MPS simulation is obtained for each value of α evaluating the maximum of the bond entropy iterating over the external field h. The bond entropy is defined by the von Neumann entropy introduced in Eq. (19). We obtain the singular values squared λ when splitting the unit cell of the infinite MPS into the bipartition. The bond entropy is one of a variety of possible entanglement measures [65,66], and it is based on the Schmidt decomposition separating a quantum system into bipartitions. The entanglement can be measured via the number of singular values, the Schmidt number, or the entropy of the singular values according to Eq. (19). For example, a Bell state has Schmidt number 2 and von Neumann entropy −2 · ( 1 2 log( 1 2 )), while an unentangled product state has Schmidt number 1 with von Neumann entropy 0. The maximum of the bond entropy is found using a grid of 21 points ranging from h [1] min = 0.9 (h [1] min = 0.0) to h [1] max = 3.9 (h [1] max = 2.0) in the ferromagnetic (antiferromagnetic) long-range Ising model. We run this evaluation for 51 different values of α ranging from α = 0 to α = 4. The critical value h [1] crit is the grid point with the maximal bond entropy. Based on this result, we refine the search on the next grid defined as where dh [i−1] is the step size on the grid. The results presented in Sec. 4 and in the following convergence studies are for h [2] crit , that is the external field at the quantum critical point between the paramagnetic and ferromagnetic, or antiferromagnetic phase. The superscript denotes the iteration refining the grid for the search and has no physical meaning past an indication for the precision of h crit .
We now discuss convergence of the iMPS results. We take the default convergence parameters of the OSMPS library [29] with the modification of an iterative bond dimension of χ = 10, 20, 40. In Fig. A2 (a)  These results call for a technical discussion of this behavior due to the large difference for many points with α 2, here in the antiferromagnetic case. This behavior is related to symmetry breaking. The degenerate ground state in the limit h → 0 is the GHZ state | ψ 0 =| ψ GHZ = |↑ · · · ↑ ± |↓ · · · ↓ √ 2 , (A.3) but without enforcing Z 2 any superposition of |↑ · · · ↑ and |↓ · · · ↓ minimizes the energy of the system. Therefore, the bond entropy can be between the bond entropy S = 0 of a product state and S = log(2) ≈ 0.69 of the GHZ ground state. We discuss the origins of this problem related to the numerical algorithm of iMPS. Tracking the singular values, we find either a distribution of pairs in the states with high entropy corresponding to the two different signs as pointed out in the limit h → 0 for the GHZ state in Eq. (A.3). The low entropy states have decaying non-degenerate singular values. The iMPS algorithm selects one of the degenerate states at non-predictable points due to the randomized initial state at the very beginning of the algorithm, as well as the entanglement truncation step, which can randomly pick a particular parity sector for the environment. This is related to the dominant eigenvector of the transfer matrix (see [20]). Note that, once one of the broken symmetry states | ↑ · · · ↑ or | ↓ · · · ↓ becomes more dominant, this dominance is "stored" for all subsequent iterations in the environment in the iMPS algorithm. Hence, symmetry breaking by random random numerical noise cannot be removed by running more iterations. The higher bond dimension and the number of iterations in the algorithm lead to smaller differences in the pairs of degenerate singular values until one of the dominant eigenvectors of the transfer matrix is kept and the other truncated. For χ = 10, the difference in the leading singular values is ≈ 0.15 decreasing to ≈ 0.05 for χ = 20. Therefore, the singular values double in magnitude since an equal counterpart was neglected. In summary we can state the following: (1) For α 1.3 symmetry breaking is not resolved. The result converges, as expected, better for faster decaying long-range interaction up to 10 −3 between the results for χ = 20 and χ = 40.
(2) For a 2 we truncate one of the symmetry broken states for χ = 40 meaning that the difference in entropy is large for any comparison to χ = 40. The results with bond dimension 20 have an error between 10 −1 and 10 −2 with regard to χ = 10. The singular values for χ = 40 are close enough to be resolved as degeneracy leading to the conclusion that the state is also well converged. (3) The region in between the algorithm alternates between resolving the symmetry broken ground state or not, leading to convergence up to 10 −4 . The same trend can be stated in the ferromagnetic case shown in Fig. A2 (b). A corresponding example for the degeneracy evolving in the singular values can be seen in Fig. A1. We emphasize that this analysis is relevant if searching for the maximal entropy on a very coarse grid or when trying to show convergence across bond dimensions. For a single bond dimension we may use the final value of the orthogonality fidelity described in Eq. (A.1). The convergence in terms of the orthogonality fidelity is discussed in Appendix C.
The finite size simulations have a two-fold purpose. On the one hand, finite size scaling provides an alternate route to obtain the thermodynamic limit of the system and is therefore a check against the iMPS results. While iMPS simulations are only limited by the bond dimension used, finite size simulations enforce the Z 2 symmetry and avoid for this reason the symmetry-breaking effects discussed in the previous section. On the other hand, we see the finite size effects affect the system, which is important as the dynamic simulations are for The key for understanding the convergence of the bond entropy is Λ, the singular values squared, shown for different bond dimensions χ. The difference in entropy is here related to the symmetry breaking appearing for χ = 40. In contrast the singular values for χ = 10 and χ = 20 appear in pairs. Considering the convergence with regards to χ by means of the bond entropy fails for this reason. This example represents a simulation for α = 4 and h = 1.14 in the ferromagnetic model. finite systems. For the finite size scaling we use the same grid method as for iMPS before with the initial grid boundaries h [1] min = 0.8 (h [1] min = 0.0) and h [1] max = 9.8 (h [1] max = 2.0) with nine (seven) grid points for each α and L in the ferromagnetic (antiferromagnetic) case. We refine the grid four times after the initial run ending with the bond dimensions χ = (256, 512) for the iteration over the OSMPS convergence parameters. The final distance between points is for both cases dh [5] ≈ 0.004. For the refinement we find h [i] crit (α, L) with the bond entropy at the middle of the system and construct the next, refined, grid i + 1 in the interval h [i] crit (α, L) ± dh [i] . The maximum number of sweeps are (4,6) for the corresponding iterations in bond dimension χ = (256, 512). The other convergence parameters, which stay constant for each set of convergence parameters iterated over, are a variance tolerance of 10 −10 and a local tolerance of 10 −10 . For power law exponents of α > 1.5, we actually reach a variance superior to 10 −7 where simulations converge better away from the critical value for the ferromagnetic Hamiltonian. In the antiferromagnetic case, longer-range interactions show worse convergence independent of the critical value of the external field. We discuss the details in Appendix B. The critical values are evaluated for system sizes 32, 64, 128, 256, and 512. By using finite size scaling, we obtain the critical value h crit of an infinite system. From Cardy [67] we know that for a finite system the distance of the critical point h crit (L) to the critical point in the thermodynamic limit h crit can be described by Knowing that the values of h crit are increasing in the ferromagnetic and antiferromagnetic case according to Fig. A3, we can use and fit the unknown parameters h crit , ν, and c via [68]. Due to the logarithmic choice of the system size, small system sizes are overrepresented in the fit. We account for this by assigning an uncertainty proportional to 1/L at each data point. The results for different system sizes and the data from the finite size scaling can be seen in Fig. A3. In the antiferromagnetic case (a), the points show the expected behavior with smaller spacing between bigger system sizes; thus the simulations are converging to the thermodynamic limit. When the interactions are not decaying for α → 0, the finite size scaling breaks down due to poorly-converging simulations. The black curve represents the L → ∞ limit achieved through the finite size scaling fit. Figure A3 (b) shows the same behavior for the ferromagnetic case. We point out that the weights for each data point are necessary in this case to avoid an intersection of the curve for L → ∞ with the finite size solution L = 512 for longer-range interactions.

Appendix A.2. Dynamic Simulations with MPS
In this section we analyze the scaling of the numerical simulation with regards to the Kibble-Zurek hypothesis. The Kibble-Zurek mechanism describes the scaling of the defect density when driving a quantum system through a quantum phase transition and was originally formulated in terms of cosmology [69], and later brought to quantum systems such as helium [70] or the nearest-neighbor quantum Ising model [26]. There, the external field is changed We point out that for the limit α → 0 and h → 0 the ground state for the ferromagnetic case is always the GHZ state, while the antiferromagnetic case has an exponential degeneracy. The latter might contribute to convergence problems when α is approaching zero. The limits α → ∞ is equal to the nearest-neighbor case and all limits h → ∞ have a ground state aligned with the external field. in order to cross the quantum critical point at h crit = 1 in the nearest-neighbor quantum Ising model. In our analysis the quantum critical point h crit depends on the actual value of the power law exponent α as presented in Sec. 2.
In order to analyze the Kibble-Zurek hypothesis in the present of long-range interactions, we first define the defect density for the ferromagnetic model. The defect density was defined in [71] using nearest-neighbor correlations, and we divide by (L − 1) instead of L to account for the open boundary condition: We now shortly outline why this works in the ferromagnetic case for the nearest-neighbor model. We recall that at the end of the quench the external field h = 0 and therefore the ground states | ψ F 0 are defined through This is equivalent to the GHZ state and we recall that this is the ground state before Z 2 symmetry breaking. The first excited state has 2(L − 1) degeneracies and is described by a single kink in the system, e.
The measurement in Eq. (A.6) yields one defect, as expected, with the spins flipping once over the chain. We emphasize that this cannot be used in any region with h = 0, except as an approximation for h < h crit . Let us turn to the defects of the ferromagnetic longrange quantum Ising model. The defects are local kinks, but the states with one kink are not all degenerate. In fact, the energy gap to the ground state is smaller for kinks closer to the boundaries. This is a manifestation of the boundary effects. Counting all kinks equally, we make an error growing in the limit α → 0. The defect density of the antiferromagnetic case in the nearest-neighbor model is defined as We again point out that the single defects have a different excitation energy depending on their position in the bulk. In order to retrieve the critical exponent scaling in the form of Eq. (23), we simulate the quenches for L = 128 for the set of τ 's containing {4, 8, 16, 32, 64, 128}. For each α we fit the coefficient c and the critical exponent µ via [68]. In order to estimate if any nonconverged simulations are included in the data, the standard deviation σ µ is calculated through the covariance matrix given by the fitting function. In Fig. A4 we show one example in the ferromagnetic model for α = 1, 2, 3, and 4. The error of the fit can be seen as an indicator for the validity of the dynamic results. Running simulations with different χ helps us to identify failing simulations versus singularities. For weaker long-range interactions with big α the fits reproduce the results while for longer-range interaction the procedure fails (α = 1). The initial state is based on the same convergence parameters as the statics discussed in Appendix A.1, except for the bond dimension which is lowered to χ = 95 or 100. For the dynamics we keep the corresponding maximal bond dimension and use the default TDVP convergence parameters in OSMPS [29] for all other parameters.

Appendix B. Convergence of Finite-Size Statics
In the appendix we briefly address the convergence of the finite size statics. We already discussed the convergence of iMPS in Appendix A and therefore it is known that the convergence can be studied in different fashions. In OSMPS, the variance of H is considered for statics of finite size systems. The user chooses a value, we take 10 −10 , and OSMPS iterates through the different bond dimension χ. We consider χ = 256 and χ = 512 of the last grid. This fifth grid has a discretization of dh ≈ 0.0044 (ferro) and dh ≈ 0.0042 (antiferro) around We start to look if the simulation converged with a variance of L·10 −10 . This information is contained in the left part of Fig. B1 for the ferromagnetic model. The coloring is as follows: green for convergence at the first convergence parameter with χ = 256, red for convergence the second convergence parameter at χ = 512, and blue if not converged at either of the previous bond dimensions. In order to estimate the variance of the blue region, we plot the actual value for the critical value h crit as a function of α and system size in the right part of Fig. B1. For the analysis used in the main part of the paper we achieve a variance less than 10 −7 .
We follow the similar scheme to study the convergence in the antiferromagnetic longrange quantum Ising Hamiltonian, presented in Fig. B2. The convergence becomes in general more difficult in the longer-range interacting region with α → 0, still reaching a variance of 10 −6 or better for the shorter-range region with α > 1.5.

Appendix C. Convergence of Infinite-Size Statics
We discussed at length that it is difficult to study convergence as a function of the bond dimension in terms of the bond entropy in Appendix A.1. But we give in the following a brief overview of the convergence of the iMPS in terms of the orthogonality fidelity defined in Eq. (A.1). In order to show convergence it is convenient to switch from the orthogonality  Figure C1 shows the convergence for the antiferromagnetic and the ferromagnetic model for the highest bond dimension χ = 40 of our convergence parameters and the second grid.
In terms of the orthogonality fidelity the simulations are well converged at an error I of 10 −8 or better for α > 1. Convergence close to the critical point is in general worse. For α < 1 interaction become too long-range to converge at the same level for these convergence settings. We remind the reader that the search interval for the ferromagnetic case shown in Fig. C1 (b) does not contain the critical value for very small α with an upper bound of h upper = 3.9.

Appendix D. Finite Size Effects in the Kibble-Zurek Scaling
Section Appendix A.2 raised the question of how finite size effects influence the value of the critical exponent µ in the dynamics. We consider the quench times τ = 4, 8, 16, 32, 64, 128 and different system sizes L = 32, 64, 128, 256 to show these finite size effects in the nearestneighbor limit. This limit corresponds to α → ∞. Table D1 presents the data. The error is solely from the fit and does not contain truncation error or methological errors from the time evolution. The thermodynamic limit is obtained via finite size scaling using the following approach: The value µ ∞ is shown in the last column of Table D1. As in the finite size scalings previous discussed, we weight the data points for large system sizes more since we have less data points in the corresponding interval. We do not show the error here, because of the limited number of data points with regards to the degrees of freedom. In conclusion, the data supports the assumption that the finite size effects are relevant and the thermodynamic limit in the nearestneighbor model approaches µ = 0.5 in both the ferromagnetic and antiferromagnetic case.  Table D1. Finite Size Effects in the Kibble-Zurek Scaling for the Quantum Ising Model.
The finite size effects lead to a growing critical exponent for smaller system sizes. In the thermodynamic limit the critical exponent of 0.5 is reached.