Lorenzo Cevolani
email: lorenzo.cevolani@institutoptique.fr

Giuseppe Carleo

Laurent Sanchez-Palencia

Spreading of correlations in exactly-solvable quantum models with

Introduction

In recent years, the study of far-from-equilibrium dynamics of correlated quantum systems has been attracting much attention [START_REF] Polkovnikov | Colloquium: Nonequilibrium dynamics of closed interacting quantum systems[END_REF][START_REF] Altman | Non equilibrium quantum dynamics in ultra-cold quantum gases[END_REF][START_REF] Eisert | Quantum many-body systems out of equilibrium[END_REF], significantly sparked by the dramatic development of experimental devices combining long coherence times, slow dynamics, and precise control of parameters. They include ultracold atoms [START_REF] Lewenstein | Ultracold atomic gases in optical lattices: Mimicking condensed matter physics and beyond[END_REF][START_REF] Bloch | Many-body physics with ultracold gases[END_REF], artificial ion crystals [START_REF] Blatt | Quantum simulations with trapped ions[END_REF], electronic circuits [START_REF] Pekola | Towards quantum thermodynamics in electronic circuits[END_REF], spin chains in organic conductors [START_REF] Giamarchi | Some experimental tests of Tomonaga-Luttinger liquids[END_REF], and quantum photonic systems [START_REF] Aspuru-Guzik | Photonic quantum simulators[END_REF]. 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 [START_REF] Kinoshita | A quantum Newton's cradle[END_REF], ballistic cone spreading of quantum correlations [START_REF] Cheneau | Light-cone-like spreading of correlations in a quantum many-body system[END_REF][START_REF] Geiger | Local relaxation and light-cone-like propagation of correlations in a trapped one-dimensional Bose gas[END_REF], thermalization effects [START_REF] Trotzky | Probing the relaxation towards equilibrium in an isolated strongly correlated one-dimensional Bose gas[END_REF][START_REF] Langen | Local emergence of thermal correlations in an isolated quantum many-body system[END_REF], nucleation of Kibble-Zurek solitons [START_REF] Lamporesi | Spontaneous creation of Kibble-Zurek solitons in a Bose-Einstein condensate[END_REF], and supercurrents in Bose superfluids [START_REF] Corman | Quench-induced supercurrents in an annular Bose gas[END_REF], 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 [START_REF] Lieb | The finite group velocity of quantum spin systems[END_REF][START_REF] Hastings | Spectral gap and exponential decay of correlations[END_REF]. 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 [START_REF] Calabrese | Time dependence of correlation functions following a quantum quench[END_REF], numerical [START_REF] Manmana | Time evolution of correlations in strongly interacting fermions after a quantum quench[END_REF][START_REF] Carleo | Light-cone effect and supersonic correlations in one-and two-dimensional bosonic superfluids[END_REF], and experimental [START_REF] Cheneau | Light-cone-like spreading of correlations in a quantum many-body system[END_REF][START_REF] Geiger | Local relaxation and light-cone-like propagation of correlations in a trapped one-dimensional Bose gas[END_REF] 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 [START_REF] Porras | Effective quantum spin systems with trapped ions[END_REF][START_REF] Deng | Effective spin quantum phases in systems of trapped ions[END_REF][START_REF] Schneider | Experimental quantum simulations of many-body physics with trapped ions[END_REF][START_REF] Islam | Onset of a quantum phase transition with a trapped ion quantum simulator[END_REF][START_REF] Kim | Quantum simulation of frustrated ising spins with trapped ions[END_REF][START_REF] Lanyon | Universal digital quantum simulation with trapped ions[END_REF][START_REF] Blatt | Quantum simulations with trapped ions[END_REF], polar molecules [START_REF] Micheli | A toolbox for lattice-spin models with polar molecules[END_REF][START_REF] Büchler | Cold Molecules: Theory, Experiment, Applications[END_REF][START_REF] Hazzard | Many-body dynamics of dipolar molecules in an optical lattice[END_REF][START_REF] Yan | Observation of dipolar spin-exchange interactions with lattice-confined polar molecules[END_REF], magnetic atoms [START_REF] Lahaye | The physics of dipolar bosonic quantum gases[END_REF], and Rydberg atoms [START_REF] Bendkowsky | Observation of ultralong-range Rydberg molecules[END_REF][START_REF] Weimer | A rydberg quantum simulator[END_REF][START_REF] Pohl | Dynamical crystallization in the dipole blockade of ultracold atoms[END_REF][START_REF] Schausz | Observation of spatially ordered structures in a two-dimensional rydberg gas[END_REF]. 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 α → ∞ [START_REF] Foss-Feig | Nearly linear light cones in long-range interacting quantum systems[END_REF][START_REF] Gong | Persistence of locality in systems with power-law interactions[END_REF][START_REF] Matsuda | Improving the Lieb-Robinson bound for long-range interactions[END_REF]. In turn, for α < D, finite-size bounds have been proposed [START_REF] Storch | Interplay of soundcone and supersonic propagation in lattice models with power law interactions[END_REF][START_REF] Kuwahara | Exponential bound on information spreading induced by quantum many-body dynamics with long-range interactions[END_REF]. 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 [START_REF] Hauke | Spread of correlations in long-range interacting quantum systems[END_REF][START_REF] Cevolani | Protected quasilocality in quantum systems with long-range interactions[END_REF] and further pointed out significant dependence to the initial state and model [START_REF] Eisert | Breakdown of quasilocality in long-range quantum lattice models[END_REF]. This is consistent with experimental observations in artificial ion chains [START_REF] Richerme | Non-local propagation of correlations in quantum systems with long-range interactions[END_REF][START_REF] Jurcevic | Quasiparticle engineering and entanglement propagation in a quantum many-body system[END_REF]. Moreover, for α > 1, the numerics showed that the propagation is significantly slower than the known bounds [START_REF] Hauke | Spread of correlations in long-range interacting quantum systems[END_REF][START_REF] Cevolani | Protected quasilocality in quantum systems with long-range interactions[END_REF][START_REF] Buyskikh | Entanglement growth and correlation spreading with variable-range interactions in spin and fermionic tunneling models[END_REF]. 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 [START_REF] Cevolani | Protected quasilocality in quantum systems with long-range interactions[END_REF], fermionic Kitaev [START_REF] Van Regemortel | Information propagation and equilibration in long-range kitaev chains[END_REF], and free-fermion [START_REF] Storch | Interplay of soundcone and supersonic propagation in lattice models with power law interactions[END_REF] 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 [START_REF] Hauke | Spread of correlations in long-range interacting quantum systems[END_REF][START_REF] Cevolani | Protected quasilocality in quantum systems with long-range interactions[END_REF], Bose [START_REF] Cevolani | Protected quasilocality in quantum systems with long-range interactions[END_REF], and Fermi [START_REF] Buyskikh | Entanglement growth and correlation spreading with variable-range interactions in spin and fermionic tunneling models[END_REF][START_REF] Van Regemortel | Information propagation and equilibration in long-range kitaev chains[END_REF] 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 [START_REF] Foss-Feig | Nearly linear light cones in long-range interacting quantum systems[END_REF][START_REF] Hastings | Locality in Quantum Systems[END_REF]. 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 [START_REF] Calabrese | Time dependence of correlation functions following a quantum quench[END_REF]. 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.

Quantum quench in quadratic Bose Hamiltonians with long-range interactions

Generic quadratic Hamiltonian

We consider a generic quadratic bosonic Hamiltonian

Ĥ = 1 2 R,R A R,R â † R âR + âR â † R + B R,R âR âR + â † R â † R , (1)
where R and R span the sites of a regular D-dimensional hypercubic lattice of unit lattice spacing. â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 [START_REF] Polkovnikov | Colloquium: Nonequilibrium dynamics of closed interacting quantum systems[END_REF]. Examples include weakly-interacting bosons and spin systems in strongly polarized states, see Refs. [START_REF] Fetter | Quantum Theory of Many-particle Systems[END_REF][START_REF] Auerbach | Interacting Electrons and Quantum Magnetism[END_REF]. 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,

B R = Uδ R + V R α (1 -δ R ) , (2)
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 as

Ĥ = 1 2 k A k â † k âk + â-k â † -k + B k â-k âk + â † k â † -k , (3)
where A k , B k , and âk are the discrete Fourier transforms of A R , B R , and âR , respectively, with the convention

f k ≡ R f R exp (ik • R), (4)
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 [START_REF] Bogoliubov | On the theory of superfluidity[END_REF],

âk = u k bk + v k b † -k , (5)
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 [ bk , b † k ] = δ k,k .
Then, provided we choose

u k = sign (A k ) 1 2 |A k | |E k | + 1 and v k = -sign (B k ) 1 2 |A k | |E k | -1 , (6)
the Hamiltonian takes the canonical form

Ĥ = E 0 + k =0 E k b † k bk , (7)
where bk and b † k represent the annihilation and creation operators of a quasi-particle of momentum k, and

E k = sign (A k ) A 2 k -B 2 k = 2sign (2h k + B k ) h k (h k + B k ) (8)
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), [START_REF] Blatt | Quantum simulations with trapped ions[END_REF], 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

|Ψ(t) = e -iH f t |Ψ 0 , (9)
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 [START_REF] Eisert | Quantum many-body systems out of equilibrium[END_REF][START_REF] Cheneau | Light-cone-like spreading of correlations in a quantum many-body system[END_REF][START_REF] Richerme | Non-local propagation of correlations in quantum systems with long-range interactions[END_REF][START_REF] Jurcevic | Quasiparticle engineering and entanglement propagation in a quantum many-body system[END_REF][START_REF] Langen | Ultracold atoms out of equilibrium[END_REF][START_REF] Geiger | Local relaxation and light-cone-like propagation of correlations in a trapped one-dimensional bose gas[END_REF][START_REF] Trotzky | Time-resolved observation and control of superexchange interactions with ultracold atoms in optical lattices[END_REF].

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

G(R, t) = Ψ(t)|(â † R + âR )(â † 0 + â0 )|Ψ(t) . (10)
For instance, this correlation function is directly connected to spin-spin correlations in the Ising model within linear spin wave theory (see Ref. [START_REF] Hauke | Spread of correlations in long-range interacting quantum systems[END_REF][START_REF] Auerbach | Interacting Electrons and Quantum Magnetism[END_REF][START_REF] Holstein | Field dependence of the intrinsic domain magnetization of a ferromagnet[END_REF] and Sec. 4.1 for details) and to the density-density correlations in the Bose-Hubbard model within mean-field theory, see Refs. [START_REF] Bogoliubov | On the theory of superfluidity[END_REF][START_REF] Natu | Dynamics of correlations in a dilute bose gas following an interaction quench[END_REF].

Turning to Fourier space and taking the thermodynamic limit it reads

G(R, t) = ˆdD k (2π) D e -ik•R Ψ 0 | â † k (t)â k (t) + â-k (t)â † -k (t) + â-k (t)â k (t) (11)
+â † 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,

âk (t) = u f k bf k (t) -v f k bf † -k (t), (12)
found from the inverse of the Bogoliubov transform [START_REF] Bloch | Many-body physics with ultracold gases[END_REF]. We then substitute the quasi-particle operator at time t by its expression bf

k (t) = exp(-iE f k t) bf k (0). ( 13
)
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

âk = u i k bi k (0) + v i k bi † -k (0) = u f k bf k (0) + v f k bf † -k (0), (14)
and

â † k = u i k bi † k (0) + v i k bi -k (0) = u f k bf † k (0) + v f k bf -k (0). ( 15
)
We then find the relation bf

k (0) = u i k u f k -v i k v f k bi k -u i k v f k -v i k u f k bi † -k . ( 16
)
Correlation spreading in bosonic quadratic Hamiltonians with long-range interactions

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 operators bi k and bi † 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 by bi 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

G c (R, t) ≡ G(R, t) -G(R, 0) (17) = 1 2 ˆdD k (2π) D F(k) e -ik•R - e i(k•R-2E f k t) + e i(k•R+2E f k t) 2
where

F(k) = A i k B f k -A f k B i k A f k + B f k E i k = h i k B i k -h f k B i k h f k + B f k E i k . ( 18
)
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 [START_REF] Lieb | The finite group velocity of quantum spin systems[END_REF][START_REF] Hastings | Spectral gap and exponential decay of correlations[END_REF]. 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 [START_REF] Calabrese | Time dependence of correlation functions following a quantum quench[END_REF]. 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,

V k = ∇ k E f k = 2 [2h k ∇ k h k + ∇ k (h k B k )] E k , (19)
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

V k = 2 E k [(2h + B k ) h + h∇ k B k ] . (20)
Conversely the parameter B k , which, according to Eq. ( 2), reads

B f k = U f + V f R =0 e ik•R |R| α , (21)
may diverge in the infrared limit, depending on the value of the exponent α. The relevant terms are

B f k -U f ∼ ˆdR ˆdΩ e ikR cos(θ) R α-D+1 ∼ k α-D and ∇ k B k ∼ k α-D-1 (22)
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.

For D + 1 < α (right column in 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. Writing

B f k = B f 0 + B f 0 |k| α-D , we find |V k | ≈ h f h f + B f 0 (α -D)B f 0 |k| D+1-α . (23)
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. Writing

B f k ≈ B f 0 /|k| D-α
, we find the energy

E f k = 2 h f B f 0 |k| D-α 2 (24)
and the velocity

|V k | = (D -α) h f B f 0 |k| D-α+2 2 . ( 25
)
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 [START_REF] Lieb | The finite group velocity of quantum spin systems[END_REF] as

G c (R, t) = g(R, t) + g ∞ (R, t) (26)
where

g ∞ (R) = 1 2 ˆdD k (2π) D e -ik•R F(k) (27)
is the asymptotic thermalization value, and

g(R, t) = - 1 2 ˆdD k (2π) D F(k) e i(k•R-2E f k t) + e i(k•R+2E f k t) 2 (28)
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

∇ k k • R ∓ 2E f k t = 0. ( 29
)
They define the separation and time dependent condition

R t = ±2∇ k E f k , (30)
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

V lc = 2Max(V k ). (31)
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 [START_REF] Lieb | The finite group velocity of quantum spin systems[END_REF]. Equation [START_REF] Büchler | Cold Molecules: Theory, Experiment, Applications[END_REF] 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,

g(R, t) λ∈Sv W (k λ ) cos k λ • R ± 2E f k λ , (32)
where the index λ spans the set of solutions of Eq. ( 30) for a fixed value of R/t, S v . The dimensiondependent quantity

W (k) = 2π t D 2 F(k) det L(k) , (33)
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) [START_REF] Cevolani | Protected quasilocality in quantum systems with long-range interactions[END_REF].

Correlation spreading in bosonic quadratic Hamiltonians with long-range interactions

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

E f k = E 0 + V 0 |k| 1-χ (34)
and the group velocity

|∇ k E f k | = (1 -χ)V 0 |k| -χ , (35)
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

G c (R, t) ∼ t R 3/2 . ( 36
)
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 [START_REF] Pohl | Dynamical crystallization in the dipole blockade of ultracold atoms[END_REF] 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 [START_REF] Hauke | Spread of correlations in long-range interacting quantum systems[END_REF][START_REF] Cevolani | Protected quasilocality in quantum systems with long-range interactions[END_REF].

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

E f k = e 0 k γ , (37)
where e 0 = 2 h f B f 0 and γ = D-α 2 . Plugging this expression into Eq. ( 17) we find

G c (R, t) ∼ ˆΩ dΩ ˆπ dkk D-1+γ e ıkR cos(θ) 1 -cos 2e 0 tk -γ , (38)
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

G c (R, t) ∼ n ı n R n n! ˆΩ dΩ cos n (θ) lim 0 + ˆπ dkk D-1+γ+n 1 + cos τ k -γ , ( 39
)
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 ˆLγ

1 π γ dqq -D+2γ+n γ [1 -cos (τ q)] = E a (-iL γ τ ) + E a (iL γ τ ) -E a (-iτ /π γ ) -E a (iτ /π γ ) 2L D+n+γ , (40)
where E a (x) is the exponential integral function of order a = D+2γ+n γ [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF]. 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

E a (-iL γ τ ) + E a (iL γ τ ) L D+n+γ ∼ sin (L γ τ ) τ 1 L 2γ+D 1 L n . ( 41
)
Plugging this expression into Eq. ( 39), we find

G c (R, t) ∼ lim L→∞ sin (L γ τ ) τ ´dΩe ı R L cos(θ) L 2γ+D . ( 42
)
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. [START_REF] Storch | Interplay of soundcone and supersonic propagation in lattice models with power law interactions[END_REF]. 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) [START_REF] Eisert | Quantum many-body systems out of equilibrium[END_REF][START_REF] Hauke | Spread of correlations in long-range interacting quantum systems[END_REF][START_REF] Buyskikh | Entanglement growth and correlation spreading with variable-range interactions in spin and fermionic tunneling models[END_REF], time-dependent variational Monte Carlo (t-VMC) [START_REF] Cevolani | Protected quasilocality in quantum systems with long-range interactions[END_REF], and discrete truncated Wigner approximation (DTWA) [START_REF] Schachenmayer | Entanglement growth in quench dynamics with variable range interactions[END_REF] and and matrix product state (MPS) [START_REF] Schachenmayer | Dynamics of correlations in two-dimensional quantum spin models with long-range interactions: a phase-space monte-carlo study[END_REF] calculations. It has also been realized experimentally using cold ion crystal chains with light-mediated long-range interactions [START_REF] Richerme | Non-local propagation of correlations in quantum systems with long-range interactions[END_REF][START_REF] Jurcevic | Quasiparticle engineering and entanglement propagation in a quantum many-body system[END_REF]. 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 [START_REF] Hauke | Spread of correlations in long-range interacting quantum systems[END_REF][START_REF] Cevolani | Protected quasilocality in quantum systems with long-range interactions[END_REF]. 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,

H = V 2 R =R σ z R σ z R |R -R | α -h R σ x R , (43)
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

H = V 2 R =R cos 2 θ σ z R σ z R + sin 2 θ σ x R σ x R -sin θ cos θ (σ x R σ z R + σ z R σ x R ) |R -R | α -h R (sin θ σ z R + cos θσ x R ) .
We then use the approximate Holstein-Primakoff transformation [52, 57]

σ z R ≈ a † R + a R and σ x R = 2n R -1 = 2a † R a R -1, (44)
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,

E cl = L D     R =R V 2|R -R | α   sin 2 θ + h cos θ   ,
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

H = E cl + V 2 R =R a † R + a R a † R + a R |R | α + 2h R a † R a R .
Up to the (irrelevant) energy E cl -2hL D , this Hamiltonian is of the quadratic form (1) with

A R,R = 2hδ R,R + V |R -R | α and B R,R = V |R -R | α . ( 45
)
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

G σσ c (R, t) = σ z R (t)σ z 0 (t) -σ z R (0)σ z 0 (0) . ( 46
)
Using the Holstein-Primakoff transformation [START_REF] Cevolani | Protected quasilocality in quantum systems with long-range interactions[END_REF], this quantity is exactly the connected correlation function defined by Eqs. ( 10) and (17).

Propagation of correlations

Figure 2 shows the space-time dynamics of the connected spin-spin correlation function G σσ c (R, t) [see Eq. ( 46)] for various values of the exponent α of the long-range exchange term and the different lattice dimensions D = 1, D = 2, and D = 3. The quench is performed by changing the value of V at fixed h in a LRTI model described by Hamiltonian [START_REF] Hauke | Spread of correlations in long-range interacting quantum systems[END_REF]. In particular we used h i = h f = 2 for almost all the quenches. The only exception is the one with α = 5/2 and D = 3, where we use h i = h f = 4 and V i = 1/2 → V f = 1/4 in order to ensure dynamical stability. For these values the dispersion relation E k is always positive and real-valued. Hence, the initial state |Ψ 0 , i.e. the ground-state of the initial Hamiltonian, is the vacuum of the magnons. For all the studied cases we checked that the condition n R 1 is fulfilled allowing the description of the LRTI model by a Hamiltonian [START_REF] Eisert | Quantum many-body systems out of equilibrium[END_REF]. These results are found by exact numerical integration of Eq. ( 17) using Eq. ( 18) for the LRTI parameters, see Eq. (4.1). In Fig. 2, the complete evolution as a function of the position R and the time t is shown in 1D, while it is plotted along the diagonals R = (R, R) in 2D and R = (R, R, R) in 3D. The results show different behaviors depending on the respective values of α and D. For D + 1 < α (right-hand-side column in 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).

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.

h i = h f = 2 and V i = 1 → V f = 1/2 in 1D

and 2D, and

h i = h f = 5 and V i = 1/2 → V f = 1/4 in 3D.
The linear system sizes are L = 2 9 in 1D, Lx = Ly = 2 8 in 2D, and Lx = Ly = Lz = 2 6 in 3D. The top panel shows the space-time dynamics of the spin-spin correlation function (color plot) together with the line R = V lc t (dashed red line)), where the light-cone velocity V lc is fitted to the boundary of the local region. The lower panel shows the comparison of the fitted light-cone velocity V lc to twice the maximum group velocity, 2 max ∂ k E f k , as computed from Eq. ( 8). Excellent agreement is found in all cases.

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. [START_REF] Hastings | Locality in Quantum Systems[END_REF], 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 [START_REF] Calabrese | Time dependence of correlation functions following a quantum quench[END_REF]. 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) = (47)
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.

h i = h f = 2, V i = 1 → V f = 1/2 and Lx = Ly = 2
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. [START_REF] Cevolani | Protected quasilocality in quantum systems with long-range interactions[END_REF]. 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 + in agreement with Ref. [START_REF] Schachenmayer | Dynamics of correlations in two-dimensional quantum spin models with long-range interactions: a phase-space monte-carlo study[END_REF], 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 [START_REF] Lieb | The finite group velocity of quantum spin systems[END_REF] 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 Correlation spreading in bosonic quadratic Hamiltonians with long-range interactions 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

τ ∼ 1 L D-α 2 . ( 48
)
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 [START_REF] Hauke | Spread of correlations in long-range interacting quantum systems[END_REF]. 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. Hence the time need to reach half the system size, L/2, scales as 1/ max(E k ) ∼ 1/L D-α 2

. G(L/2; L/2; τ * ) Moreover, the analytic approach used in Sec. 3.3 also provides the scaling of the amplitude of the correlation function at t = τ * . It yields

G σσ c (L/2, τ * ) ∝ τ * L 2γ+D = 1 L 3D-α 2 (49)
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 [START_REF] Carleo | Light-cone effect and supersonic correlations in one-and two-dimensional bosonic superfluids[END_REF]. 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 [START_REF] Hauke | Spread of correlations in long-range interacting quantum systems[END_REF][START_REF] Cevolani | Protected quasilocality in quantum systems with long-range interactions[END_REF]. 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 [START_REF] Hauke | Spread of correlations in long-range interacting quantum systems[END_REF] and t-VMC [START_REF] Cevolani | Protected quasilocality in quantum systems with long-range interactions[END_REF] 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 [START_REF] Hauke | Spread of correlations in long-range interacting quantum systems[END_REF][START_REF] Cevolani | Protected quasilocality in quantum systems with long-range interactions[END_REF]. So far causality breaking has been observed experimentally in one-dimensional ion chains of moderate sizes [START_REF] Richerme | Non-local propagation of correlations in quantum systems with long-range interactions[END_REF][START_REF] Jurcevic | Quasiparticle engineering and entanglement propagation in a quantum many-body system[END_REF]. 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 [START_REF] Kim | Quantum simulation of frustrated ising spins with trapped ions[END_REF][START_REF] Lanyon | Universal digital quantum simulation with trapped ions[END_REF][START_REF] Richerme | Non-local propagation of correlations in quantum systems with long-range interactions[END_REF][START_REF] Jurcevic | Quasiparticle engineering and entanglement propagation in a quantum many-body system[END_REF][START_REF] Britton | Engineered two-dimensional ising interactions in a trapped-ion quantum simulator with hundreds of spins[END_REF][START_REF] Porras | Quantum manipulation of trapped ions in two dimensional coulomb crystals[END_REF], polar molecules [START_REF] Yan | Observation of dipolar spin-exchange interactions with lattice-confined polar molecules[END_REF][START_REF] Pupillo | Condensed Matter Physics with Cold Polar Molecules[END_REF], magnetic atoms [START_REF] Menotti | Dipolar interaction in ultra-cold atomic gases[END_REF][START_REF] Lu | Strongly dipolar bose-einstein condensate of dysprosium[END_REF][START_REF] Aikawa | Bose-einstein condensation of erbium[END_REF][START_REF] Paz | Nonequilibrium quantum magnetism in a dipolar lattice gas[END_REF], Rydberg atoms [START_REF] Schausz | Observation of spatially ordered structures in a two-dimensional rydberg gas[END_REF][START_REF] Macrì | Rydberg dressing of atoms in optical lattices[END_REF][START_REF] Schwarzkopf | Imaging spatial correlations of rydberg excitations in cold atom clouds[END_REF][START_REF] Schwarzkopf | Spatial correlations between rydberg atoms in an optical dipole trap[END_REF][START_REF] Anderson | Dephasing of resonant energy transfer in a cold rydberg gas[END_REF][START_REF] Nipper | Highly resolved measurements of stark-tuned förster resonances between rydberg atoms[END_REF][START_REF] Butscher | Atom-molecule coherence for ultralong-range rydberg dimers[END_REF], and alkaline Earth atoms [START_REF] Olmos | Long-range interacting many-body systems with alkaline-earth-metal atoms[END_REF][START_REF] Rey | Probing many-body interactions in an optical lattice clock[END_REF][START_REF] Martin | A quantum many-body spin system in an optical lattice clock[END_REF][START_REF] Lemke | p-wave cold collisions in an optical lattice clock[END_REF][START_REF] Swallows | Suppression of collisional shifts in a strongly interacting lattice clock[END_REF]. It is expected that the analysis in terms of diverging quasi-particle energies and group velocities is generic to all systems. However, the boundaries between the local, quasi-local, and non-local regimes can be affected by long-range terms in the free component of the quadratic Hamiltonian. Moreover special care should be taken on the relative weights of local and nonlocal contributions to the correlation functions. For instance, it has been shown that causality is protected irrespective to the strength of the long-range interactions in the extended Bose-Hubbard model in dimension D = 1 [START_REF] Cevolani | Protected quasilocality in quantum systems with long-range interactions[END_REF]. It would be interesting to study the behavior of the correlation function of the same model in dimensions higher than one.

Figure 1 :

 1 Figure 1: Energy (top panels) and modulus of the group velocity V f k = |V f k | (bottom panels) for a two-dimensional (D = 2) long-range system described by the dispersion relation (8) with h k = h, for various values of the exponent α. For instance, this dispersion relation emerges in the long-range transverse Ising (LRTI) model discussed in Sec. 4.1.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.

2 Figure 2 :Figure 3 :

 223 Figure 2: Space-time evolution of the spin-spin correlation function following quantum quenches in the LRTI model in various dimensions D (rows) and for different values of the exponent α (columns).For almost all cases, the quenches are from V i = 1/2 to V f = 1 for a fixed magnetic field h = 2. The only exception is for the 3D case with α = 5/2 (left, bottom) where we used the quench V i = 1/4 → V f = 1/2 and h i = h f = 4 in order to avoid dynamical instabilities. The linear system sizes are L = 2 12 in 1D, Lx = Ly = 2 9 in 2D, and Lx = Ly = Lz = 2 6 in 3D, with periodic boundary conditions. Distances are measured in units of the lattice constant and times in units of the inverse magnetic field.

Figure 4 :

 4 Figure 4: Top panel: t as function of R different values of in log-log scale (points) with the fitted function (continuous lines). It is possible to see the algebraic dependence of t from R and the agreement between the fit and the numerical data. Bottom panel: Function β ( ) as extracted from fits to t (R) = R β for different values of . It is possible to see that as becomes smaller β ( ) approaches a constant, black line. The points in red and blue correspond to the parameters obtained by the fit of the data of the same colors in the top panel. The data used for this analysis comes from a quench in a D = 2 model with α = 2.3 andh i = h f = 2, V i = 1 → V f = 1/2and Lx = Ly = 2 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.

Figure 5 :

 5 Figure 5: Values of the fit parameter β as a function of α for systems in dimensions D = 1 and D = 2. The data are obtained analyzing the time evolution of correlation in systems of length L = 2 13 = 8192 and L = 2 6 = 64 for D = 1 and Lx = Ly = 10 for D = 2. The inset in the left figure presents the dependence of the parameter β on the system size L for the case of α = 3/2 and D = 1.

Figure 6 :

 6 Figure 6: Activation time (upper panels) and amplitude (lower panels) of the spin-spin correlation function computed at R = L/2 in the non local regime (α < D) for a 1D (left column) and 2D (right column) systems of different sizes. Note the log-log plot scales. The net decrease of the time of the first maximum for different values of α and L is the clear signature of locality breaking. The numerical data (points) are in good agreement with the analytical predictions (straight lines). The slopes of the straight lines is fixed by Eq. (48) and their intercepts have been found fitting the numerical data .

Acknowledgments

We acknowledge the discussion with P. Hauke. This research was supported by the European Research Council (FP7/2007-2013 Grant Agreement No. 256294), Marie Curie IEF (FP7/2007-2013 Grant Agreement No. 327143), and FET-Proactive QUIC (H2020 Grant No. 641122). It was performed using HPC resources from GENCI-CCRT/CINES (Grant No. c2015056853). Use of the computing facility cluster GMPCS of the LUMAT Federation (FR LUMAT 2764) is also acknowledged.

In this appendix we detail the calculations for the scaling of the correlations function G (R, t) in the case D = 1 and α = 3/2 outlined in Sec. 3.2. We start the computation from the general point of view, ie for generic D and D < α < D + 1 where the correlation function take form

) where E 0 and V 0 are finite constants coming from the expansion of the dispersion relation around k ≈ 0. The correlation horizon is expected to be determined by the contributions with largest velocities, that is around the divergent point k = 0, we can write the correlation function [START_REF] Lieb | The finite group velocity of quantum spin systems[END_REF] as

We then focus on the first integral and write it as a power series in t ˆπ 0 dk cos (kR

3) This new integral can be computed for every value of R. We find

where 1 F 2 the hypergeometric function [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF]. For large values of R, we use the asymptotic limit of the latter, which yields ˆπ 0 dk cos (kR)

We can evaluate simply the summation over n of the first term

In the limit of large R and t the last term will go to zero leaving the correlation function unaffected.

Inserting now B 1 n in Eq. (A.3), we need to compute the sum over n and this is possible analytically just for χ = 1/2, which corresponds to D = 1 and α = 3/2. For these values of the parameters we find

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

The integral can be again expressed as a hypergeometric function for every value of n ˆπ 0 cos (kR)

taking now the asymptotic value of this function in the large R limit we get ˆπ 0 cos (kR)

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

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 [START_REF] Lieb | The finite group velocity of quantum spin systems[END_REF] thus scales as

(A. [START_REF] Langen | Local emergence of thermal correlations in an isolated quantum many-body system[END_REF] as discussed in the main text.


Antoine Amarilli

Silviu Maniu

Mikaël Monet

Challenges for Efficient Query Evaluation on Structured Probabilistic Data

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

To have an accurate description of the real world, it is often necessary to associate probabilities to our observations. For instance, experimental and scientific data may be inherently uncertain, because, e.g., of imperfect sensor precision, harmful interferences, or incorrect modelling. Even when crisp data can be obtained, it can still be the case that we do not trust who retrieved it or how it came to us. The notion of probabilistic databases has been introduced to capture this uncertainty, reason over it, and query it: these databases are augmented with probability information to describe how uncertain each data item is. Given a probabilistic database D and a query q, the probabilistic query evaluation problem (PQE) asks for the probability that the query q holds on D. Unfortunately, even on the simplest probabilistic database models, PQE is generally intractable [START_REF] Dalvi | Efficient query evaluation on probabilistic databases[END_REF].

One possibility to work around this intractability is to use approximate approaches, such as Monte Carlo sampling on the data instances. A different direction was recently explored in [START_REF] Amarilli | Provenance circuits for trees and treelike instances[END_REF], namely, restricting the kind of input instances that we allow, in what we call the structural approach. It is shown in [START_REF] Amarilli | Provenance circuits for trees and treelike instances[END_REF] that the data complexity of PQE is linear if the instances have bounded treewidth, i.e., they can be structurally decomposed in a tree-like structure where each node must contain at most k elements, for a fixed parameter k. Moreover, in [START_REF] Amarilli | Tractable lineages on treelike instances: Limits and extensions[END_REF], it is shown that bounding the instance treewidth is necessary to ensure the tractability of PQE, because some queries are hard on any unbounded-treewidth instance family (under some conditions). Hence, bounded-treewidth methods seem to be the right way to make PQE tractable by the structural approach.

These theoretical works, however, left open the question of practical applicability: many challenges must still be addressed to implement a practical system using these techniques. First, obtaining an optimal decomposition of an arbitrary instance is NP-hard [START_REF] Arnborg | Complexity of finding embeddings in a k-tree[END_REF]. Second, the complexity is only polynomial in the data, with the query and parameter being fixed; this hides a constant which can be exponential in the width k and non-elementary in the query q. Third, we do not Sec. 5

Fig. 1: Overview of the structural approach for PQE know which real datasets can indeed be decomposed, at least partially, with a small k. This paper thus presents our vision of a database management system based on the structural approach, and gives an overview of the research directions, both theoretical and practical, which we intend to address to this end.

Probabilistic Query Evaluation: A Structural Approach

We first review our structural approach [START_REF] Amarilli | Provenance circuits for trees and treelike instances[END_REF] for probabilistic query evaluation (PQE). The approach is illustrated in Figure 1.

The approach applies to tuple-independent (TID) instances (but generalizes to more expressive models [START_REF] Amarilli | Leveraging the Structure of Uncertain Data[END_REF]). Formally, a TID instance I is a relational database D where each tuple t ∈ D has some probability p t . The TID I represents a probability distribution over the subinstances D ⊆ D (subsets of facts): following the independence assumption, the probability of D is t∈D p t × t∈D\D (1 -p t ).

We study the probabilistic query evaluation (PQE) problem: given a Boolean query q and TID instance I, determine the probability that q holds on I, i.e., the total probability of the subinstances of I that satisfy q. We refer to the combined complexity of PQE when I and q are given as input; we refer to data complexity when I is the input and q is fixed.

The first step of the structural approach (Section 3) is to translate the query q to a formalism that can be efficiently evaluated. In the approach of [START_REF] Amarilli | Provenance circuits for trees and treelike instances[END_REF], following [START_REF] Courcelle | The monadic second-order logic of graphs. I. Recognizable sets of finite graphs[END_REF], the query is compiled to a tree automaton, i.e., a finite-state automaton over trees [START_REF] Comon | Tree automata: Techniques and applications[END_REF]. The approach works for expressive queries written in monadic second-order logic, which covers in particular first-order logic and (unions of) conjunctive queries. This translation of the query is independent from the instance, so does not affect data complexity; however, it depends on a parameter k of the instance, to be defined soon. Intuitively, the automaton represents an algorithm to evaluate the query on suitable instances.

The second step (Section 4) applies to the instance I, and computes a structural decomposition of it. In [START_REF] Amarilli | Provenance circuits for trees and treelike instances[END_REF], we compute a tree decomposition [START_REF] Bodlaender | Treewidth computations I. Upper bounds[END_REF][START_REF] Bodlaender | Treewidth computations II. Lower bounds[END_REF], equivalent to junction trees in graphical models [START_REF] Lauritzen | Local computations with probabilities on graphical structures and their application to expert systems[END_REF], and then a tree encoding over a finite alphabet: the results of [START_REF] Amarilli | Tractable lineages on treelike instances: Limits and extensions[END_REF] show that tree decompositions are essentially the only possible way to make PQE tractable. The parameter k measures how well I could be decomposed: in our case, k is the treewidth, measuring how close I is to a tree. By treelike instances, we mean instances whose treewidth is bounded by a constant.

The third step (Section 5) is to compute a lineage of the query q on the instance I, i.e., compute an object that represents concisely the subinstances of I that satisfy q. This object can be used for PQE, as what we want to compute is precisely the total probability of this set of subinstances. Specifically, we compute a Boolean lineage circuit of the tree automaton for the query over the tree encoding of the instance, according to the construction of [START_REF] Amarilli | Provenance circuits for trees and treelike instances[END_REF]. This step is purely symbolic and does not perform any numerical probability computation.

The fourth and last step is to evaluate efficiently the probability of the query from this lineage representation, by computing the probability that the circuit is true. This task cannot be performed efficiently on arbitrary Boolean circuits, but it is feasible in our context, for two independent reasons [START_REF] Amarilli | Provenance circuits for trees and treelike instances[END_REF][START_REF] Amarilli | Tractable lineages on treelike instances: Limits and extensions[END_REF]. First, this circuit can also be tree decomposed, which allows us to apply a messagepassing algorithm [START_REF] Lauritzen | Local computations with probabilities on graphical structures and their application to expert systems[END_REF] for efficient probability computation. Second, in the case where we made the query automaton deterministic [START_REF] Comon | Tree automata: Techniques and applications[END_REF], the circuit is actually a d-DNNF [START_REF] Darwiche | On the tractable counting of theory models and its application to truth maintenance and belief revision[END_REF], for which probabilistic query evaluation is tractable.

Efficient Compilation to Expressive Automata

Compiling the query to an automaton following the structural approach of [START_REF] Amarilli | Provenance circuits for trees and treelike instances[END_REF], by applying [START_REF] Courcelle | The monadic second-order logic of graphs. I. Recognizable sets of finite graphs[END_REF], is generally non-elementary in the query. This section presents our main ideas to address this problem: we intend to restrict to tractable query fragments, and to use more expressive automata targets to compile the query more efficiently. These challenges are not specific to PQE; the next section presents the lineage computation tasks, which are specific to PQE.

Efficient Compilation. Of course, we cannot hope to compile all Monadic Second Order logic (MSO) to automata efficiently, or even all conjunctive queries (CQs). Indeed, efficient compilation to automata implies that non-probabilistic query evaluation is also tractable in combined complexity on treelike instances; however, CQs are already hard to evaluate in this sense (even on fixed instances). Hence, we can only hope to compile restricted query languages efficiently.

Many fragments are known from earlier work to enjoy efficient combined query evaluation. In the database context, for instance, acyclic CQs can be evaluated in polynomial combined complexity [START_REF] Yannakakis | Algorithms for acyclic database schemes[END_REF]. This generalizes to first-order logical sentences that can be written with at most k variables, i.e., FO k [START_REF] Grädel | On logics with two variables[END_REF]. However, it also generalizes to the guarded fragment (GF) [START_REF] Andréka | Modal languages and bounded fragments of predicate logic[END_REF], whose combined complexity is also PTIME, and where better bounds can be derived if we know the instance treewidth [START_REF] Berwanger | Games and model checking for guarded logics[END_REF]. The tractability of GF, however, does not capture other interesting query classes: reachability queries, and more generally two-way regular path queries (2RPQs) and variants thereof [START_REF] Baeza | Querying graph databases[END_REF], as well as Monadic Datalog as in [START_REF] Gottlob | Monadic Datalog over finite structures of bounded treewidth[END_REF].

Our first task would thus be to develop an expressive query language that captures GF, 2RPQs, and Monadic Datalog. Ideally this fragment should be parameterized, i.e., all CQs or all FO queries q could be expressible in the fragment, up to increasing some parameter k q , with the compilation being PTIME for fixed k q but intractable in k q . We would then develop an efficient algorithm to compile such queries to automata that check them on bounded-treewidth instances, for fixed values of the query parameter k q and of the treewidth. Our ongoing work in this direction investigates very recent extensions of GF with negation and fixpoints [START_REF] Benedikt | A step up in expressiveness of decidable fixpoint logics[END_REF][START_REF] Benedikt | The complexity of boundedness for guarded logics[END_REF], for which compilation to automata was studied as a tool for logical satisfiability. We believe that these results, suitably extended and adapted to query evaluation, can yield to bounded-treewidth automaton compilation methods that covers the query classes that we mentioned.

Expressive Automata Targets. The efficient compilation of queries to automata is made easier by allowing more expressive automaton classes as the target language. In [START_REF] Amarilli | Provenance circuits for trees and treelike instances[END_REF], we used bottom-up tree automata, which process the tree decomposition of the instance from the leaves to the root. Our idea is to move to more expressive representations, namely, two-way alternating automata [START_REF] Comon | Tree automata: Techniques and applications[END_REF]. These automata can navigate through the tree in every direction (including already visited parts), and thus can be more concise. The notion of alternation allows automata to change states based on complex Boolean formulae on the neighboring states, which also helps for concision. Indeed, the expressive languages of [START_REF] Benedikt | A step up in expressiveness of decidable fixpoint logics[END_REF] are compiled to two-way alternating parity automata, which further use a parity acceptance condition on infinite runs, to evaluate fixpoints.

To make automaton compilation more efficient, another idea is to compile queries to automata with a concise implicit representation. In particular, we can use automata with a structured state space: the states are tuples of Boolean values, and the transition function can be written concisely for each coordinate of the tuple as a function of the tuples of child states. It may be possible to capture the tractability of query evaluation for 2RPQs via automaton methods, structuring the state space to memorize separately the regular sublanguages of paths between node pairs.

Obtaining Tree Decompositions

Estimating Treewidth. As we have mentioned, computing the treewidth of an instance directly is an NP-hard problem. Hence, a practical system using the structural approach must compute tree decompositions more efficiently, even if this limits us to non-optimal decompositions. We intend to experiment with two main kinds of methods to obtain tree decompositions efficiently: separator-based algorithms, which recursively divide the instance based on various heuristics; and elimination ordering algorithms, where the nodes in the graph are ordered using some measure and eliminated one by one from the graph [START_REF] Bodlaender | Treewidth computations I. Upper bounds[END_REF]. To estimate the quality of our decompositions, we can also estimate lower bounds on the instance treewidth: for instance via graph degeneracy or average degree [START_REF] Bodlaender | Treewidth computations II. Lower bounds[END_REF].

Query-Specific Decompositions. In some cases, knowledge about the query can help us to obtain better tree decompositions of the instance. A trivial situation is when we know that the query is only on a subset of the database relations: we can then ignore the others when decomposing. More subtly, if we know that specific joins are not made by the query, then we may be able to rewrite the instance accordingly, and lower the treewidth. For instance, if no R-and S-atoms share a variable in the query, then the instance {R(a, b), S(b, c)} can be rewritten to {R(a, b), S(b , c)}, which may lower the treewidth by disconnecting elements. We do not understand this process yet in the general case, but we believe that a theory of lineage-preserving instance rewritings for a given query (or query class) can be developed, using the notion of instance unfoldings introduced in [START_REF] Amarilli | Tractable lineages on treelike instances: Limits and extensions[END_REF].

Tractable Lineage Targets

Once we have compiled the query to an automaton and decomposed the instance to a encoding, our goal is to compute a lineage representation of the automaton on the encoding, namely, a representation of the subinstances where the query holds, which we will build as a Boolean circuit. We can then use this to perform PQE, by computing the probability of the query as that of the lineage. In so doing, we need to rely on the fact that the lineage is in a class of circuits for which probability can be efficiently computed.

To this end, a first step towards a practical system is to adapt the methods of [START_REF] Amarilli | Provenance circuits for trees and treelike instances[END_REF] to the expressive automaton classes that are needed for efficient query compilation. We believe that this is possible, but with a twist: because two-way automata can navigate a tree in every direction, they may go back from where they came, thus resulting in cyclic runs. Therefore, it seems that the natural lineages that we would obtain for alternating two-way automata are cyclic Boolean circuits, which we call Boolean cycluits. A semantics for such circuits would need to be defined based on the semantics of automaton runs and reachable states: we believe that the evaluation could follow least fixed-point semantics, and be performed in linear time.

Second, we would need to perform efficient probability computation on these cycluits. One first way to address this would be to eliminate cycles and transform them to tractable classes of Boolean circuits (e.g., d-DNNFs), which we believe can be done assuming bounds on the treewidth of the cycluits. Alternatively, we can apply message passing methods directly on the cycluits [START_REF] Lauritzen | Local computations with probabilities on graphical structures and their application to expert systems[END_REF]; or we can try to rewrite the automaton to produce acyclic circuits or even d-DNNFs directly. All these methods would be generally intractable in the query, which is unsurprising: indeed, PQE is often intractable even for languages with tractable combined complexity, and efficient compilation to automata. It would be interesting, however, to identify islands of tractability; and, in intractable cases, to benchmark the previously mentioned approaches and see which ones perform best in practice.

Another important direction for a practical system is to be able to evaluate queries on instances where facts are not independent, i.e., go beyond the TID formalism. For instance, facts could be present or absent according to a complex lineage, like the cc-instances of [START_REF] Amarilli | Leveraging the Structure of Uncertain Data[END_REF]. In this context, new methods can be efficient. For instance, if the number of probabilistic events is small, performing Shannon expansions on some well-chosen events may make large parts of the instance deterministic, making the query easier to evaluate on these parts.

We now review possible approaches and directions to implement and evaluate the structural approach for PQE on real-world datasets.

Results on Treelike Instances. In [START_REF] Monet | Probabilistic evaluation of expressive queries on bounded-treewidth instances[END_REF], the structural approach has been compared with one of the existing probabilistic data management systems, namely MayBMS [START_REF] Huang | MayBMS: a probabilistic database management system[END_REF]. The instances considered have been randomly generated to have low treewidth [START_REF] Benedikt | A step up in expressiveness of decidable fixpoint logics[END_REF]. The results show that an implementation of the structural approach can perform query evaluation faster than the exact methods of MayBMS, in cases where there are many matches and many correlations between them. Indeed, MayBMS does not take advantage of the fact that the instances are treelike. However, in this work, the queries were compiled to automata by hand rather than automatically, and there was no study of practical datasets.

Practical Datasets and Partial Decompositions. A first question is to extend this study to practical datasets, and to investigate whether such datasets have low treewidth, or whether we can use approximate decompositions or reasonably low treewidth. Our preliminary results suggest that some datasets have high treewidth, but others, in particular transportation networks, have treewidth much smaller than their size. For instance, the OpenStreetMaps graph of Paris has over 4 million nodes and 5 million edges, but we estimated its treewidth to be 521. We do not know yet of a theoretical reason explaining why transportation networks generally exhibit this property.

However, this bound is still too large to be practical. One way to work around this problem is thus to compute a partial decomposition [START_REF] Wei | TEDI: Efficient shortest path query answering on graphs[END_REF] of the instance, i.e., a tree decomposition of a part of the instance whose width is at most k, with k fixed. This results in a structure formed of a forest of instances with treewidth k, called the tentacles, that interface with a core, i.e., the remaining facts, whose treewidth is too high and that cannot be decomposed. Our preliminary experiments have shown that, for some transportation networks, a partial decomposition for k = 10 results in a core instance whose size is about 10% of the original instance.

This decrease in the size of the core, in turn, can potentially have an immediate effect in the processing of queries. Preliminary results [START_REF] Maniu | ProbTree: A query-efficient representation of probabilistic graphs[END_REF] have shown that using partial decompositions of fixed treewidth for probabilistic reachability queries, in conjunction with sampling in the core graph, can make query processing up to 5 times faster.

Tentacle Summarization. An important problem when computing probabilities on partial decompositions is the interface between the tentacles and the core, i.e., we must find a way to summarize the tentacles in the core when applying sampling to the core. As the tentacles are treelike, we can efficiently compute probabilities and lineages in them: the goal of summarization is to eliminate the tentacles and replace them by summary facts that are added to the core. In the case of simple queries, such as reachability queries, the summary facts can have the same semantics as in the original instance, but this does not seem to generalize to arbitrary queries: it may even be the case that some queries cannot be rewritten to the summary facts while remaining in the same language.

Having summarized the tentacles, we may also answer queries approximately via sampling, using the (exact) tentacle summaries added to the core: as the instance is now smaller, this process can be performed faster.