Spreading of correlations in exactly-solvable quantum models with long-range interactions in arbitrary dimensions

We study the out-of-equilibrium dynamics induced by quantum quenches in quadratic Hamiltonians featuring both short- and long-range interactions. The spreading of correlations in the presence of algebraic decaying interactions, $1/R^\alpha$, is studied for lattice Bose models in arbitrary dimension $D$. These models are exactly solvable and provide useful insight in the universal description of more complex systems as well as comparisons to the known universal upper bounds for the spreading of correlations. Using analytical calculations of the dominant terms and full numerical integration of all quasi-particle contributions, we identify three distinct dynamical regimes. For strong decay of nteractions, $\alpha>D+1$, we find a causal regime, qualitatively similar to what previously found for short-range interactions. This regime is characterized by ballistic (linear cone) spreading of the correlations with a cone velocity equal to twice the maximum group velocity of the quasi-particles. For weak decay of interactions, $\alpha<D$, we find instantaneous activation of correlations at arbitrary distance. This signals the breaking of causality, which can be associated with the divergence of the quasi-particle energy spectrum. Finite-size scaling of the activation time precisely confirms this interpretation. For intermediate decay of interactions, $D<\alpha<D+1$, we find a sub-ballistic, algebraic (bent cone) spreading and determine the corresponding exponent as a function of $\alpha$. These outcomes generalize existing results for one-dimensional systems to arbitrary dimension. We precisely relate the three regimes to the first- and second-order divergences of the quasi-particle energy spectrum for any dimension.


Introduction
In recent years, the study of far-from-equilibrium dynamics of correlated quantum systems has been attracting much attention [1][2][3], significantly sparked by the dramatic development of experimental devices combining long coherence times, slow dynamics, and precise control of parameters. They include ultracold atoms [4,5], artificial ion crystals [6], electronic circuits [7], spin chains in organic conductors [8], and quantum photonic systems [9]. In ultracold-atom systems for instance, major assets are the possibility to engineer out-of-equilibrium initial states and dynamically change some microscopic parameter(s) of the system. The study of the system dynamics after such a quantum quench now makes it possible to address a variety of open basic questions with unprecedented accuracy. So far, it paved the way to observation of undamped oscillations in integrable onedimensional systems [10], ballistic cone spreading of quantum correlations [11,12], thermalization effects [13,14], nucleation of Kibble-Zurek solitons [15], and supercurrents in Bose superfluids [16], for instance.
Universal properties of the time evolution of local observables following a quantum quench mainly relies on so-called Lieb-Robinson bounds. For lattice systems with short-range interactions, the correlation function of any set of two local observables can be activated only after some finite time t , which increases linearly with the distance R between the supports of the two observables [17,18]. This provides an upper bound to the spreading velocity of correlations. It corresponds to a cone in space-time representation, which defines a causal region. In many cases, the known bounds provide a fair account of the actual spreading of correlations for shortrange interactions. Ballistic (linear cone) behavior has now been found in several analytical [19], numerical [20,21], and experimental [11,12] works.
The extension of the Lieb-Robinson bounds to quantum systems with long-range interactions is a major challenge, with possible applications to a variety of systems, including artificial ion crystals [22][23][24][25][26][27][28], polar molecules [29][30][31][32], magnetic atoms [33], and Rydberg atoms [34][35][36][37]. An important step forward in the understanding of the dynamics of lattice systems with two-body long-range interactions of the form 1/R α was made with the identification of logarithmic Lieb-Robinson-like bounds, t ∼ log (R), for strong-enough decay of interactions, α > D. It was further shown that for α > 2D, the bounds can be made more stringent in the form of a polynomial-shaped horizon, t ∼ R β , where β smoothly converges to β → 1 (linear cone) for α → ∞ [38][39][40]. In turn, for α < D, finite-size bounds have been proposed [41,42]. However, no bound is known in the thermodynamic limit, which suggests possible instantaneous activation of correlations at arbitrary distance, and correspondingly, the breaking of causality. Numerical work confirmed the breaking of causality in one-dimensional lattice spin models for α < 1 [43,44] and further pointed out significant dependence to the initial state and model [45]. This is consistent with experimental observations in artificial ion chains [46,47]. Moreover, for α > 1, the numerics showed that the propagation is significantly slower than the known bounds [43,44,48]. More precisely, it was found to be sub-ballistic for 1 < α < 2 and ballistic for α > 2. Finally, non universal behavior was found in certain systems. For instance, in the extended one-dimensional Bose-Hubbard [44], fermionic Kitaev [49], and free-fermion [41] chains with long-range interactions, clear ballistic spreading was found irrespective to the interaction exponent α, which corresponds to efficient dynamical protection of causality in these systems.
In view of this rich behavior, analyzing exactly-solvable systems is thus of utmost importance to determine the precise dynamics of quantum correlations beyond mathematically-exact bounds, which are not guaranteed to be saturated. In this respect, quadratic Hamiltonians play a central role. For instance, quadratic approximations have been studied for one-dimensional spin [43,44], Bose [44], and Fermi [48,49] systems. In the present work, we consider quadratic Bose systems in arbitrary lattice dimension D, hence generalizing previous results to dimensions higher than one. We develop the general theory of correlation dynamics for Bose systems undergoing an instantaneous quantum quench between two quadratic Hamiltonians with both short-and longrange interactions of the form 1/R α . We provide the equations for the time evolution of a generic correlation function, which can be easily generalized to more complicated cases. Then we study the first-and second-order divergences of the energy spectrum as a function of α and D, and precisely relate them to the dynamical behavior of the correlations by computing analytically the dominant contributions. For strong decay of interactions, α > D + 1, the group velocity of the quasi-particle excitations is bounded, which yields a linear conic causal region. This behavior is similar to that found for short-range interactions and corresponds to a dynamics significantly slower than the known bounds [38,50]. For weak decay of interactions, α < D, the energy spectrum diverges in the infrared limit. It provides a vanishing characteristic time, independent of the distance R, for the activation of correlations. The latter can be associated with instantaneous propagation of correlations and the breakdown of causality. This is compatible with the absence of known finite bound in this regime. Finite-size scaling of the typical times precisely confirms this behavior. For intermediate decay of interactions, D < α < D + 1, we find a bent-cone causal region determined by a sub-ballistic algebraic bound, t ∼ R β , where β is some function of the exponent α and the dimension D. This again corresponds to a dynamics that is significantly slower than the known bounds. Furthermore, we study the specific long-range transverse Ising model in dimensions D = 1, 2, and 3 in the (quadratic) spin-wave approximation. We study the full space-time dynamics of the spin-spin correlations for various values of α. Taking into account the contributions of all quasiparticles, we confirm the three regimes. We then characterize each regime in detail. For α < D, we perform finite-size scaling of the correlation function, which confirms our analytical predictions for both the bound and the amplitude of the correlations at the propagation front, and the breaking of causality. For α > D + 1, we find a clear linear cone. We determine the associated velocity and find excellent agreement with the excepted value of twice the maximum group velocity [19]. For D < α < D +1, we find a clear algebraic bound, t ∼ R β for all tested cases and extract numerically the exponent β(α) in dimensions D = 1 and 2. Finally we study the shape of the correlation front in dimension D > 1 and discuss its symmetries.

Generic quadratic Hamiltonian
We consider a generic quadratic bosonic Hamiltonian where R and R span the sites of a regular D-dimensional hypercubic lattice of unit lattice spacing. a R andâ † R are, respectively, the annihilation and creation operators at site R, with the usual bosonic commutation relations [â R ,â † R ] = δ R,R , and the coefficients A R,R and B R,R are coupling amplitudes, containing both short-and long-range terms. A variety of systems can be described by the Hamiltonian (1). Examples include weakly-interacting bosons and spin systems in strongly polarized states, see Refs. [51,52]. Without loss of generality, we write A R = 2h R + B R . For simplicity, we assume h R is short range, while the long-range character of interactions is entirely included in the coefficient B R . More precisely, we assume that it contains a contact interaction term and an algebraic long-range interacting term, where U is the on-site contact interaction strength, V is the long-range interaction strength, and α is some non negative constant. Generalization to the case where h R also contains long-range interactions is straightforward. Assuming translation invariance and parity symmetry, the coefficients A R,R and B R,R only depend on the Cartesian inter-site distance R = |R − R |. This condition allows us to write Hamiltonian (1) in momentum space aŝ where A k , B k , andâ k are the discrete Fourier transforms of A R , B R , andâ R , respectively, with the convention for any field f R . The annihilation and creation operatorsâ k andâ † k fulfill the bosonic commutation rule [â k ,â † k ] = δ k,k and, due to parity symmetry, the coefficients A k and B k are real-valued. Hamiltonian (3) can now be diagonalized using the standard Bogoliubov transformation [53], where the functions u k and v k can be assumed to be real-valued without loss of generality and to fulfill condition u 2 k − v 2 k = 1 to ensure the commutation relation [b k ,b † k ] = δ k,k . Then, provided we choose the Hamiltonian takes the canonical form whereb k andb † k represent the annihilation and creation operators of a quasi-particle of momentum k, and is the quasi-particle dispersion relation. The quantity E 0 is the zero-point energy, i.e. the energy of the vacuum of quasi-particles. Dynamical stability requires that the quasi-particle energy E k is real-valued, i.e. h k (h k + B k ) ≥ 0. Equations (5), (6), and (8) provide the complete information to determine any equilibrium and out-of-equilibrium properties of the system.

Quantum quench and correlation function
We focus our attention on the out-of-equilibrium dynamical properties of the system induced by a quantum quench. This protocol consists in preparing the system in some initial state |Ψ 0 at time t = 0 and let it evolve under the action of some final Hamiltonian H f . For instance, |Ψ 0 may be the ground state of another initial Hamiltonian H i . Here we assume that H i and H f are both generic quadratic bosonic Hamiltonians of the form of Eq. (1) and the quench amounts to an abrupt change of the amplitudes A R and B R from the values A i R and B i R to the values A f R and B f R . Assuming that the quench H i → H f takes place on a time scale shorter than any characteristic dynamical time, the time evolution of the system for t > 0 is determined by the equation where we set = 1. Quantum quenches constitute a controlled protocol to study out-of-equilibrium dynamics of correlated quantum systems and are now experimentally realised in ultracold-atom systems [3,11,46,47,[54][55][56].
The post-quench dynamical properties of the system can be studied via the correlation function of some local observableΣ R . Here we consider the simplest observable that can be constructed from the local annihilation and creation operators, that isΣ R =â R +â † R . The corresponding correlation function is then For instance, this correlation function is directly connected to spin-spin correlations in the Ising model within linear spin wave theory (see Ref. [43,52,57] and Sec. 4.1 for details) and to the density-density correlations in the Bose-Hubbard model within mean-field theory, see Refs. [53,58]. Turning to Fourier space and taking the thermodynamic limit it reads +â † k (t)â † −k (t) |Ψ 0 in the Heisenberg picture. In order to compute explicitly the correlation function G(R, t), we first substitute the particle annihilation and creation operators by their expressions in terms of the quasi-particle ones associated to the final Hamiltonian, found from the inverse of the Bogoliubov transform (5). We then substitute the quasi-particle operator at time t by its expression We finally use the two Bogoliubov transformations (5) associated to the initial and final Hamiltonians respectively. At the time of the quench, t = 0, they yield We then find the relation This expression allows us to write the correlation function G(R, t) as a function of the position R and of the time t, and the initial quasi-particle operatorsb i k andb i † k . We then calculate the quantum average over the initial state |Ψ 0 , which we assume to be the ground state of the initial Hamiltonian H i , defined byb i k |Ψ 0 = 0 for any k. After some straightforward algebra we find that the correlation function can be evaluated in the thermodynamic limit. It reads The first right-hand-side term in Eq. (17) is the asymptotic thermalization value while the second right-hand-side term is the time-dependent part. The latter contains the information on the spreading of correlations in the system we are interested in.

Relevant divergences due to long-range interactions
Before analyzing the dynamical behavior of the correlation function G c (R, t) found above, it is worth discussing the divergences that appear in the various terms of Eq. (17) due to long-range interactions. This is motivated by the known dynamical behavior of short-range systems. For shortrange interacting quantum systems the propagation of correlations following a quantum quench exhibits a light-cone structure in its space-time dynamics [17,18]. It shows up in the form of a linearly increasing horizon, which can be interpreted as being generated by the contribution of the fastest quasi-particles defined by the final Hamiltonian H f . For that class of Hamiltonians, the velocity defined by the horizon is then expected to be twice the maximum group velocity of the quasi-particles [19]. This is expected to hold in general, including for long-range systems, whenever the post-quench Hamiltonian has excitations with well-defined, finite group velocities. In contrast, sufficiently long-range interactions can make the group velocity diverge and the quasi-particle picture breaks down [43-45, 49, 59] possibly yielding non ballistic propagation of correlations. Divergences in the energy spectrum, not only in the group velocity, may further affect the dynamics of correlations. When E f k diverges, at a finite k value, the space coordinate becomes irrelevant and the associated characteristic time τ ∼ 1/E f k vanishes. This may yield instantaneous activation of correlations at arbitrary distance and, consequently the breaking of locality. Note that this scenario is not incompatible with the quasi-particle picture, since divergence of the energy E f k at finite k implies divergence of the group velocity V k = ∇ k E f k and the break-down of the quasi-particle picture.
It is thus expected that the relevant divergences are those of the quasi-particle spectrum, Eq. (8), and of the quasi-particle group velocity, of the post-quench Hamiltonian. We further assume that the parameter h k is regular and it admits a linear development h k ≈ h + h · k in the infrared limit. We then obtain the limit expression Conversely the parameter B k , which, according to Eq. (2), reads may diverge in the infrared limit, depending on the value of the exponent α. The relevant terms are in the thermodynamic limit where Ω is the D dimensional solid angle and θ is the azimuthal angle with respect to k. Typical behaviors of the energy and group velocity for various values of α are shown in Fig. 1. Fig. 1), both B k and ∇ k B k converge to a finite value in the infrared limit. Hence, both the energy and the group velocity are bounded for any value of the momentum k. Note that the maximum group velocity is not necessary at k = 0 as for instance in the example shown in the figure.
For D < α < D + 1 (central column in Fig. 1), B k is finite but ∇ k B k diverges in the infrared limit. Hence, the group velocity diverges, giving rise to infinitely fast modes, while the energy is finite with a cusp around the origin.
For α < D (left column in Fig. 1), both B k and ∇ k B k diverge in the infrared limit. Hence, both the energy and the group velocity diverge.
and the velocity As shown in the following section, those various behaviors play a central role in the qualitative space-time behavior of the correlation function. For α < D (left panels), both the energy and the group velocity diverge in k = 0. For D < α < D + 1 (central panels) the energy is finite but shows a cusp around k = 0, which corresponds to a divergent group velocity around the same point. For D + 1 < α, both the energy and the group velocity are finite and well behaved. Note that the absolute maximum of the group velocity is located close to but not exactly at the origin k = 0.

Horizon shape as a function of α
We now show how the properties of the post-quench excitation spectrum E f k and the weight F(k) determine the locality horizon and its breaking. It results from the discussion of Sec. 2.3 that we have to distinguish three cases with different divergence properties.

Local Regime (D + 1 < α)
Consider first the case where both the energy and the group velocity are bounded in the whole Brillouin zone for the quadratic Hamiltonian described in Sec. 2.1. This occurs for algebraically decaying interactions of the type Eq. (2) with α < D + 1.
To study the evolution of the correlation function, it is worth separating the static and timedependent components, and rewrite the correlation function (17) as where is the asymptotic thermalization value, and is the time-dependent part. The latter contains the relevant time dependence of the correlation function given by the post-quench dynamics. This contribution may be interpreted as the spreading of two counter-propagating beams of quasi-particles, which are represented by the two oscillating functions e i(k·R∓2E f k t) . Using the stationary-phase approximation, the main contribution to Eq. (28) is given by the stationary points of the two phase terms, determined by the two equations They define the separation and time dependent condition where the ± sign represents the two directions of the beams. This procedure can be interpreted as selecting the contribution to the correlation function (28) of the modes with a velocity equal to R/t. Since the group velocity V k = ∇ k E f k is bounded for α < D + 1, it has a maximum value V M . Then Eq. (30) has solutions only for |R|/t ≤ V M . This defines a ballistic (linear) horizon, that is a "light cone", in the |R| − t plane. Its slope gives the "light-cone" velocity V lc , defined by The presence of a ballistic horizon in the out-of-equilibrium dynamics is thus directly connected to the presence of a finite absolute maximum of the group velocity [17]. Equation (30) can also predict what happens for points outside the "light-cone". If |R|/t exceeds the maximum value 2V M , then Eq. (30) has no solution. In this case the integration over the oscillating functions has no stationary point and the correlation functions is suppressed. More precisely for |R|/t < 2V M the contribution to the time-dependent part of the correlation function of the modes with parameter v = R/t is given by the stationary-phase-approximation expression in generic dimension, where the index λ spans the set of solutions of Eq. (30) for a fixed value of R/t, S v . The dimensiondependent quantity where L D is the determinant of the Hessian matrix of the final dispersion relation E f k , represents the weight associated to each contributing pair of modes. In practice, some of the modes with the velocity R/t may be insignificant if they have an extremely small weight W(k λ ) compared to the other modes with the same velocity. This circumstance, however, does not affect the horizon, as long as at least one mode has a significant weight. In the opposite case the effective spreading of correlations may be slower than the expected bound, Eq. (31) [44].

Quasi-local regime (D < α < D + 1)
Let us now turn to the case where the energy is finite but the group velocity diverges due to a cusp in the energy spectrum around k = 0. It follows from Eq. (8) and the discussion of Sec. 2.3 that the dispersion relation of the post-quench Hamiltonian may be written and the group velocity where χ = D + 1 − α. An analytic approach is extremely difficult because of the complexity of the computations. In Appendix A we detail the computation for the case D = 1 and α = 3/2, where an analytic result can be computed. The result involves the explicit computation of the contribution to the correlation function coming from the modes around k ≈ 0, where the group velocity diverges.The correlation function scales algebraically in the large R and t region scales as Hence, the correlation horizon is algebraic, t ∼ R β . While this result is found here for a specific case, we show numerically in Sec. 4.2.2 that it is true for all tested values of the dimension D and of the exponent α. The case analyzed in this section gives β = α, as we will see in Sec. 4.2.2 this cannot be generalized to other values of α where in general, we numerically find β = α. Note that the scaling (36) is slower than ballistic. This is surprising because it may be expected that a divergent group velocity would allow faster-than-ballistic scaling. This idea is also consistent with extended Lieb-Robinson bounds, which are faster-than-ballistic. Our analysis shows that interference effects between the contributing divergent modes strongly affect the correlation front and the known bounds are not saturated. Note that similar behavior was found in 1D spin chains using many-body numerical approaches [43,44].

Non-local regime (α < D)
Consider finally the case where the quasi-particle energy spectrum (1) diverges. As discussed in Sec. 2.2 this is the case for α < D, owing to the divergence of the Fourier transform of the potential (2). The dispersion relation around k = 0 takes the form where e 0 = 2 h f B f 0 and γ = D−α 2 . Plugging this expression into Eq. (17) we find the factor k γ comes from the contribution of the weight F ∼ 1/E i k . Since the integral is dominated by the low-k components, the upper bound π of the integral is irrelevant. The lower bound k = holds for finite-size systems of linear length L and scales as ∼ 1/L. Hence the limit → 0 is equivalent to the thermodynamic limit L → ∞. We proceed by expanding the previous expression in powers of R and find where we use the dimensionless time τ = 2e 0 t. We can then integrate this expression term by term using the transformation k → q = k −γ and find where E a (x) is the exponential integral function of order a = D+2γ+n γ [60]. In the last expression, the last two terms are bounded and the limit L → ∞ can be taken without any problem after the summation. We thus focus on the first two terms, which contain the diverging energy contributions that will affect locality. In the large R limit, we find Plugging this expression into Eq. (39), we find The last equation shows that an algebraic divergence in the quasi-particles spectrum can lead to a signal which appears on a time scale 1/L γ and it goes to zero in the thermodynamic limit. Note that this time scale is directly connected to the divergence in the energy spectrum with the same exponent γ. This parameter depends on the specific model and interaction decaying. In our case γ = D−α 2 while the same results applies to free-fermion chain with γ = D 2 − α. The latter is consistent with the scaling found in Ref. [41]. In this regime, the function sin (L γ τ ) /τ gives rise to a contribution of the type δ (τ ) to the correlation function at any distance. The same expression can be used to obtain the scaling of the value of the correlation function itself as G (R, t) ∼ 1/L γ+D . Moreover, these expressions show that the dominant contributions to the correlation function carry spherical symmetry despite the underlying lattice geometry. This will be important for our discussion of the correlation front in Sec 4.3. In the next section we will check all the analytic prediction made in this and in the last section.

Application to the long-range transverse Ising model
In order to compare the analytic predictions of Sec. 3 to exact correlation functions for a physically relevant quadratic Hamiltonian, we now consider a specific case, namely the long-range transverse Ising (LRTI) model. The one-dimensional version of the latter has been studied previously in the presence or absence of a transverse field using time-dependent density matrix renormalization group (t-DMRG) [3,43,48], time-dependent variational Monte Carlo (t-VMC) [44], and discrete truncated Wigner approximation (DTWA) [61] and and matrix product state (MPS) [62] calculations. It has also been realized experimentally using cold ion crystal chains with light-mediated long-range interactions [46,47]. Spin wave analysis within quadratic approximation has been shown to correctly predict the dynamical regimes and shows good agreement with full numerical calculations in the polarized regime [43,44]. Here we extend this analysis to an arbitrary lattice dimension D.

Spin-wave analysis of the long-range transverse Ising model
Let us first briefly recall the quadratic approximation for the LRTI model, where σ j R with j ∈ {x, y, z} are the local Pauli matrices and |R − R | is the Cartesian distance between the two sites R and R on the D-dimensional hypercubic lattice. In order to write Hamiltonian (43) into the quadratic form (1) we use linear spin wave theory (LSWT). We first rotate the reference axes around the free axis y by an arbitrary angle θ. In the rotated frame, the new spin operators read σ x R = cos θ σ x R − sin θ σ z R , σ y R = σ y R , and σ z R = sin θ σ x R + cos θ σ z R , and the Hamiltonian We then use the approximate Holstein-Primakoff transformation [52,57] σ valid for small bosonic occupation number n R 1, and expand the Hamiltonian in the form H = n≥0 H n where every H n contains exactly n Holstein-Primakoff particle operators among a R , a R , a † R , and a † R . The zeroth-order term is the classical energy, where L D is the total number of lattice sites. The rotation angle is chosen to minimize the classical energy, which yields θ = π for antiferromagnetic exchange, V > 0. The Hamiltonian computed in the configuration with θ = π then reads Up to the (irrelevant) energy E cl − 2hL D , this Hamiltonian is of the quadratic form (1) with The LRTI model is thus of the form discussed in Sec. 2.1 with a parameter h k = h that does not depend on the momentum k. In particular, the dispersion relation E k = 2 h (h + B k ) is a monotonous decreasing function of the modulus of the momentum, k = |k|.
In the following, we consider the time-dependent dynamics of the connected spin-spin correlation function Using the Holstein-Primakoff transformation (44), this quantity is exactly the connected correlation function defined by Eqs. (10) and (17).  Fig. 2), the results show clear evidence of a strong form of locality, namely ballistic cone spreading. While the correlations are significant for t > R/V lc , where V lc is some velocity, they are instead strongly suppressed for t < R/V lc . For D < α < D + 1, we still find evidence of locality with correlations appearing for t > F (R), where F is some finite-valued function. This behavior is clear in 1D and 2D while, in 3D, finite lattice size effects hardly permits to determine the function F (see details below). For α < D, the numerical data is compatible with locality breakdown and instantaneous activation of the correlations, irrespective of the distance. Still a very thin band with vanishing correlations is visible at short times. It is due to finite-size effects and their scaling actually confirms locality breakdown (see details below).

Propagation of correlations
This behavior is qualitatively consistent with the previous analysis for the different regimes. In the following, we discuss the three identified regimes and provide a quantitative comparison between the analytic predictions and the numerical data.   4.2.1. Local Regime (D + 1 < α) -Let us start with the local regime, which corresponds to fast decay of the long-range exchange term, α > D + 1. In this regime, the stationary-phase analysis predicts ballistic spreading of the correlations with a velocity equal to twice the maximum group velocity (see Sec. 3.1). This is confirmed by the complete numerical data. The upper panels of Fig. 3 show the space-time dynamics of the correlation function, similarly to Fig. 2, but in contour plot format and with a strong color contrast. It shows a clear ballistic (light-cone-like) behavior of the correlation front in all dimensions. Fitting a linear function, R = R 0 + V lc t, to the correlation front, we find the light-cone velocity V lc . The results are compared to the predicted value of twice the maximum group velocity of the final Hamiltonian in the lower panels of Fig. 3. We find an excellent agreement for all the studied cases within the error bars. The width of the error bars reflects the fact that the leaks outside the light-cone are algebraically decaying, Ref. [50], and this makes more complicated to define the exact position of the correlation front. The good quantitative agreement between the numerics and the prediction of Eq. (31) confirms that the correlation front is mainly determined by the propagation of counter-propagating quasi-particles with the highest velocities, whenever they exist, as predicted by the Cardy-Calabrese scenario [19].

4.2.2.
Quasi-local regime (D < α < D + 1) -In Sec. 3.2 we have demonstrated that the correlation front for α = 3/2 in D = 1 scales as t/R 3/2 but, due to the complexity of the calculation, it has not been possible to extend this result to other cases. It is then important to study the behavior of the correlation horizon for different values of α and different dimensions D. We impose a threshold and then we find the first value of time t when the correlation function reaches for every value of R,Ḡ (t , R) = In particular, we consider time-averaged correlation functions,Ḡ (t, R) = 1 t´t 0 dτ G σσ c (R, t), in order to minimize the effects of undesirable small time oscillations.  11 . Such system size is necessary to get a good fit in the large R region, where the algebraic regime is supposed to be found.
In the top panel of Fig. 4 the values of t as function of R for a D = 2 system are shown for different values of in log-log scale. From these plot, it is clear that there is an algebraic dependence between these two quantities in the large R regime, as suggested by the analytic result for a specific case in Sec. 3.2. We can then interpolate these points with a generic algebraic dependence of the type t (R) = t 0 + m * R β for every values of . The limit → 0 + will give us the correct and independent scaling of the horizon lim →0 β( ). This limit can be found in the bottom panel of Fig. 4 where the values of the fitted parameter β are plotted as function of for α = 2.3, and it is possible to see how these results corrects the ones obtained in Ref. [44]. In Fig. 5 the values of β as functions of α are plotted for D = 1 and D = 2 systems of different sizes. It is possible to see that β → 1 as α → D + 1, in agreement with Ref. [62], which means that there is a continuous transition between the non ballistic, D < α < D + 1, and the ballistic, α > D + 1, regimes. On the other side, the transition at α = D between the non -local and the non ballistic regime, is discontinuous. From our data, it is possible to extrapolate the two limits. For D = 1, we find β = 1.52 ± 0.02 for α → 1 and β = 1.01 ± 0.08 for α → 2.For D = 2, we fnd β = 1.56 ± 0.3 for α → 2 and β = 1.1 ± 0.5 and α → 3. This can be explained directly from the expression (17) and from the divergences studied in Sec. 3. In the region α < D the dispersion relation is explicitly divergent, and this leads to the non-local regime studied in Sec. 3.3. For all the values α > D the dispersion relation itself is not divergent and depends continuously on α, which means that in this region the function β has to be continuous too. This motivates the discontinuity of the function β for in α = D and its continuity in α = D + 1. This last point can be explained naively saying that approaching α = D +1 the divergence in the derivative of the dispersion relation disappears, leaving the spectrum without divergences. We now discuss finite-size effects, which are important as we will see. In Fig. 5, we show a comparison of the values of the parameter β for two 1D systems of different sizes, namely L = 2 13 = 8192 and L = 2 6 = 64. In spite of coresponding to system sizes that differ by more than two orders of magnitude, the results are quite close. In particular, they yield α → 1 and α → 2 limits that are consistents within error bars. Nevertheless, the results for the largest system are systematically above those found for the smallest system. In order to get more insight on finite-size effects, we have studied the behavior of β versus the system size for α = 3/2 and D = 1. The results shown on the Inset of Fig. 5 show a systematic increase up to the largest system size we are able to compute. It shows that very large systems are necessary to reach the thermodynamic limit. However, the value of β we find for L = 2 13 is β 1.34, which is in fair agreement with the analytic prediction, β = 1.5 within 10%.

4.2.
3. Non local regime (α < D) -We finally discuss the non local regime, which corresponds to very weak algebraic decay of the exchange interaction, with α < D. The breaking of locality is apparent in the plots of Fig. 2 (left column). According to the discussion of Sec. 3.3, it may be attributed to the infrared divergence of the energy spectrum, which corresponds to a vanishing typical activation time of correlations at arbitrary distance in the thermodynamic limit. In order to corroborate the estimate of Sec. 3.3, we may take advantage of finite-size effects. In this case the minimal time scale is provided by the inverse of the maximum energy scale, which corresponds to the momentum k ∼ 1/L. This can be used to obtain an estimate of the system-size dependent activation time The latter determines the arrival time of the first maximum of the correlation function for large distances. In Fig. 6(a) and (b) we plot the arrival time τ * of the first maximum of the spinspin correlation function at a distance equal to half the system size, R = L/2, as a function of L in 1D and 2D. Excellent agreement between the numerical data (points) and the predicted scaling (48) (dashed lines) is found for various values of α < D in both 1D and 2D. These results confirm the predictions of Sec. 3.3. Note that the same scaling can be found from the quasi-particle approach [43]. For power-law spectra as considered here, the group velocity scales as V k E k /k, whose maximum is found for k ∼ 1/L. Moreover, the analytic approach used in Sec. 3.3 also provides the scaling of the amplitude of the correlation function at t = τ * . It yields Figures 6(c) and (d) compare numerical data (points) to the analytic prediction above (dashed lines) for the amplitude of the correlation function at R = L/2 and τ * . Good agreement is found in both 1D and 2D. The outcome further confirms the analytic predictions of Sec. 3.3. Note that this result is a direct consequence of the interference between the fastest modes and cannot be found by the simplest independent quasi-particle approach.  In this section we finally discuss the shape of the propagation front during the time evolution of the correlation function. For simplicity, we consider the two-dimensional case but the extension to higher dimensions D is straightforward. In Fig. 7 the correlation function G σσ c (R, t), Eq. (46), for a quench in a D = 2 LRTI system is plotted as function of the position R at various times t and for different values of the exponent α. In the non causal regime, α < D, the correlation function is significantly different from zero for every values of R at any time t. Conversely, in the causal and quasi-causal regimes, for α > D, correlations take a finite time to be activated, which increases as a function of the distance R. For α > D + 1, a sharp edge is visible in the correlations and it evolves ballistically in time. In contrast, for D < α < D + 1, the correlation front has a different scaling which is the signature of the quasi-locality. This is consistent with the discussion of Sec. 4.2.
Let us now focus on the correlation pattern. For α < D the correlation function is spherically symmetric for large values of R while in the region close to the origin this symmetry is no longer present. This is in perfect agreement with Eq. (3.3) which predicts a spherical symmetry of the correlation function in the large R region. For α > D + 1 there is a well-defined correlation front that spreads in the system and its symmetry is spherical despite the presence of the lattice. The symmetry of the front is due to the fact that the maximum group velocity is located very close to k = 0, where the spectrum is spherically symmetric (see Fig. 1). The inner structure of the correlation function is due to the other two local maxima, which are not in the infrared region and whose contribution to the correlation function is not spherically symmetric. This contrasts with the behavior observed for the short-range Bose-Hubbard model, where the maximum group velocity is located at finite k and gives rise to a non spherical correlation front in 2D [21]. For the quasi-local regime D < α < D + 1 we can use the same arguments used for the other two regimes. The divergence in the velocity located at k = 0 is not sufficient to destroy completely locality as discussed in Sec. 3.3 and a sort of locality, called quasi-locality, appears. Still, as for the other two regimes, the modes that dominate close to the horizon are located in the infrared region, where the spherical symmetry of the long-range potential dominates. This determines the spherical symmetry of the correlation function in the large R region. These considerations can be extended straightforwardly to any dimension higher than one because they only rely on the analysis of the symmetries of the energy spectrum and in particular around the point where is located the maximum group velocity that dominates the evolution of the correlation horizon in all regimes.

Conclusions
In this work we have studied the space-time spreading of correlations in a bosonic quadratic Hamiltonian with long-range interactions in hypercubic lattices of arbitrary dimension D. We have assumed that the interaction term decays algebraically with the exponent α. The dynamics is induced by an instantaneous quench of the Hamiltonian parameters at the initial time. We have shown that the spreading of correlations is determined by the first-and second-order divergence properties of the energy spectrum of the well-defined quasi-particles, i.e. the divergences of the energy and the group velocity, hence generalizing previous results available in dimension D = 1 [43,44]. We have introduced a generic expression for the space-time evolution of the correlation function. We have identified three distinct regimes in the spreading of correlations. In the case where the quasi-particle energy and group velocity is finite (α > D+1), the dynamics shows a strong form of causality, characterized by a ballistic spreading of correlations. The propagation velocity, so-called light-cone velocity, is determined by the propagation of quasi-particles of opposite and maximum velocity, and is thus equal to twice the maximum group velocity. This behavior is equivalent to what happens in short-range interacting systems. In the case where the quasi-particle energy is finite but the group velocity diverges (D < α < D + 1), the space-time behavior of the correlation function instead results from the interference of the quasi-particle contributions with high velocities. This yields a non-ballistic correlation front. The latter was found to be algebraic, t ∼ R β , and sub-ballistic, β > 1, in all studied dimensions D and exponents α. This is consistent with and extends previous numerical calculations using t-DMRG [43] and t-VMC [44] calculations performed in dimension D = 1. In the case where the quasi-particle energy diverges, the activation of correlations is instantaneous, hence leading to complete breaking of causality. This can be attributed to a vanishing activation time in the thermodynamic limit. We have provided an analytical formula for the finite-size scaling of the activation time and correlation function, which confirms the breaking of causality in this system in any dimension.
Our analytic predictions are supported by the complete calculation of the space-time dynamics of the correlation function for the bosonic quadratic Hamiltonian corresponding to the linear spin wave approximation of the long-range transverse Ising model in dimensions D = 1, D = 2, and D = 3, as well as by many-body numerical approaches in dimension D = 1 [43,44]. So far causality breaking has been observed experimentally in one-dimensional ion chains of moderate sizes [46,47]. Our results pave the way to the experimental observation of causality and its breaking in dimensions higher than one. Several atomic, molecular, and optical systems exhibit long-range interactions, which can be controlled. They include artificial crystal ions [26,27,46,47,63,64], polar molecules [32,65], magnetic atoms [66][67][68][69], Rydberg atoms [37,[70][71][72][73][74][75], and alkaline Earth atoms [76][77][78][79][80]. It is expected that the analysis in terms of diverging quasi-particle energies and This function scales as t/R 3/2 multiplied by a smooth oscillating function.
The second term in Eq. (A.2) can be studied along the same lines. First we write it as a power series in t sin (2E 0 t)ˆπ As we demonstrated in Eq. (A.7), the summation of A 2 n over n goes to zero as 1/R and it will not affect the correlation function in the regime defined by large R and t. We can plug B 2 n in Eq. (A.9) and sum over n. As before, it is possible to perform these computations analytically in the case D = 1 and χ = 1/2 and it gives − sin (2E 0 t) ∞ n=0 cos π 4 (2n + 1) Γ 1 + 1 2 (2n + 1) (2V 0 t) 2n+1 R 1+ 1 2 (2n+1) (2n + 1)!
Again, this term scales as t/R 3/2 and the oscillating functions do not affect this main behavior. For D = 1 and α = 3/2, the correlation function (17) thus scales as G c (R, t) ∼ t R 3/2 . (A.14) as discussed in the main text.