Discrete truncated Wigner approach to dynamical phase transitions in Ising models after a quantum quench

By means of the discrete truncated Wigner approximation we study dynamical phase transitions arising in the steady state of transverse-field Ising models after a quantum quench. Starting from a fully polarized ferromagnetic initial condition these transitions separate a phase with nonvanishing magnetization along the ordering direction from a symmetric phase upon increasing the transverse field. We consider two paradigmatic cases, a one-dimensional long-range model with power-law interactions $\propto 1/r^{\alpha}$ decaying algebraically as a function of distance $r$ and a two-dimensional system with short-range nearest-neighbour interactions. In the former case we identify dynamical phase transitions for $\alpha \lesssim 2$ and we extract the critical exponents from a data collapse of the steady state magnetization for up to 1200 lattice sites. We find identical exponents for $\alpha \lesssim 1$, suggesting that the dynamical transitions in this regime fall into the same universality class as the nonergodic mean-field limit. The two-dimensional Ising model is believed to be thermalizing, which we also confirm using exact diagonalization for small system sizes. Thus, the dynamical transition is expected to correspond to the thermal phase transition, which is consistent with our data upon comparing to equilibrium quantum Monte-Carlo simulations. We further test the accuracy of the discrete truncated Wigner approximation by comparing against numerically exact methods such as exact diagonalization, tensor network as well as artificial neural network states and we find good quantitative agreement on the accessible time scales.


I. INTRODUCTION
Recent impressive developments underline the rich phase structures that can be generated by forcing isolated quantum matter out of equilibrium. Some examples of these phenomena are the emergence of exotic phases, loss of adiabaticity across critical points in the context of the Kibble-Zurek mechanism and non-equilibrium phase transitions. These are some of the multiple aspects currently at the centre of an intense theoretical and experimental activity, as summarized in the reviews [1][2][3][4][5][6][7].
A paradigmatic protocol to drive a many-body system out of equilibrium, routinely used in experiments and intensively studied theoretically, is a quantum quench. After initializing the system in a state, that can be thought of as a ground state of a given initial Hamiltonian, it is let evolving after an abrupt change of a Hamiltonian parameter. The long-time steady states after such quantum quenches can feature symmetry-broken phases and singular behaviour at the transition towards the disordered phase. These Dynamical Phase Transitions (DPTs) [8][9][10] may be understood as transitions in the micro-canonical ensemble in case the many-body system thermalizes, driven by shifting the system's energy across the symmetry-restoration threshold. In non-ergodic systems, however, long-time steady states can be realized which cannot be described in terms of the conventional thermodynamic ensembles. As a particular consequence, such systems allow the generation of phases and phase transitions with properties that cannot be realised in any equilibrium context [11].
In this work we focus on DPTs realized in spin-1/2 Ising models in transverse fields. We consider the case of long-range interacting models, which have recently attracted a lot of attention [8,10,[12][13][14][15][16][17][18][19][20] and constitute a paradigmatic class of non-ergodic systems capable of generating non-equilibrium steady states as reported both theoretically [8][9][10]21] and experimentally [22][23][24]. It was shown [8,10] that starting from an initial fully polarized state along the ordering direction, the asymptotic state of these systems can undergo a transition from an ordered phase at small fields to a disordered one when the field exceeds a critical value. While inherently of nonequilibrium character, the resulting phases can be characterized by means of the conventional Landau paradigm via local order parameters. Still, the understanding of the nature of the transition between the ordered and disordered phases has remained limited. In particular, it is unclear to which extent these DPTs follow the general paradigm of continuous equilibrium transitions such as to whether they can be categorized in terms of universality classes and therefore whether the concepts of universality and scaling extend to this non-equilibrium dynamical regime. We remark that here we completely neglect the analysis of singular behaviours in the (infinite-size) time dynamics, which is another aspect of DPTs [5] with some connection with the symmetry-breaking behaviour [10].
In this work we show that the DPTs after a quantum quenches in transverse-field Ising chains with power-law decaying interactions (∼ r −α ) can feature scale invari-ance. We find evidence that the critical exponents of the DPT are universal over a large range of interaction exponents α. Via finite-size scaling of the time-averaged longitudinal magnetization we identify the critical value of the field h c of the DPT and, in particular, determine the scaling exponents of the transition. By studying the decay in time of the longitudinal magnetization we are able to put bounds to the values of α above which the ordered phase disappears. We can confirm the existence of two phases as long as α 2. The critical field drops rapidly to zero for values of α larger or equal than 2. This indicates that only the trivial phase survives in this regime of α consistent with previous works [10].
For α = 0 the dynamics can be solved via an effective mean field description, which becomes exact in the thermodynamic limit. For α ≥ 0 we compute the quantum real-time evolution by means of the Discrete Truncated Wigner Approximation (DTWA) [25]. It has already been reported that DTWA compares well with other methods for long-range models [26,27] and, as we are going to show, works very well also for our problem, giving a very good comparison with the results of a recent numerical study using tensor network methods [10]. The DTWA has the advantage that it allows us to access large sizes with moderate computational resources polynomially scaling in the system size. Consequently, we can perform finite-size scaling also in long-range systems where it is crucial to reach large system sizes in order to tell the difference from the infinite-range (α = 0) case.
When analyzing scale invariance at the DPT, we find that the DTWA gives rise to scaling exponents identical to the mean-field ones at α = 0. For finite α, at the mean-field level, the exponents are of course independent on the range of the interaction. This is different for the DTWA, it compares well with exact methods, as emphasized above, and it is able to capture correlations. Therefore, in principle, it can give reliable scaling exponents. We computed the dependence on α of the scaling exponents of the magnetization, and observed a significant deviation from the mean-field values at α ∼ 1. As discussed in the relevant sections, in this regime of α DTWA is not able to achieve accurate precision for a reliable scaling. It clearly indicates however when the deviations from mean field occur.
The favourable scaling with the DTWA with the number of sites allows to tackle the study the DPT also in higher dimensions, a problem never touched so far in the literature. As long as the spins are interacting via longrange exchange couplings we do not expect significant dependence on the dimensionality. This is why we decided to study a two-dimensional system with nearestneighbour coupling. Also in this case we expect a DPT. Here however, the critical behaviour should clearly deviate from the mean-field case. In this case we can only compare to exact diagonalization at small system sizes to test the quality of the DTWA approach. As we will show in the second part of the paper, we are able to detect the existence of the DPT through an analysis of the magnetization and of the Binder cumulant. We perform a comparison with finite-temperature quantum Monte Carlo results and we see that this DPT corresponds to the thermal transition. We show that this result is physically sound because the model is quantum chaotic. We find additional support for this conclusion by using exact diagonalization and showing that the level-spacing statistics is Wigner-Dyson.
The paper is organized as follows. In Section II we introduce the model and also define the order parameter for the phase transition. In Section III we discuss the DTWA theory in detail and we show how to apply it to our model. In Section III B we compare the DTWA approach for this model with known results both in the infinite-range interaction case -where exact diagonalization is possible also for large sizes -and long-range interaction where the TDVP method is used. For the range of parameters we are interested in, we find that the comparison is very good. In Section IV A we perform the finitesize scaling analysis for the one-dimensional long-range case. We first consider the case α = 0 where we compare with the exact diagonalization results and find that the comparison is very good. Then we move to analyze the case α = 0 and see that the transition exists only when α 2. The results for the short-range two-dimensional models are discussed in Section IV B. Finally, Section V is devoted to the conclusions and further perspectives. The appendices contain additional details of the numerical analysis.

II. THE MODEL
As anticipated in the introduction, we will study a system of N interacting spins governed by the Hamiltonian where theσ x,z j are the Pauli matrices of the spin located in the j−th site, h is an external transverse field and J ij is the exchange coupling between the spins. We will consider two cases (assuming to express the energies in units of the exchange coupling): • A long-range interacting spin exchange in one dimension. We assume periodic boundary conditions defining the distance between two sites as r ij = min[|i − j|, N − |i − j|]. The Kac factor ensures that the Hamiltonian is extensive.
• A short-range interacting spin on a d-dimensional cubic lattice where the exchange coupling is different from zero (δ l,m is the Kronecker-delta) only if i and j are nearest-neighbours (nn). We will assume periodic boundary conditions and consider the cases of d = 1 and d = 2 (a square lattice of size L, N = L 2 ).
The system is initialized in the state fully polarized along x, We then perform a quantum quench with the dynamics governed by the Hamiltonian of Eq. (1). We are interested in the evolution of the total x (longitudinal) magnetization which is given by and the order parameter for the DPT is the long-time average of this magnetization.
In the short-range two-dimensional case, we will also analyze the Binder cumulant in the long-time limit, defined as where we defined m (l) x i l |ψ(t) . The Binder cumulant is a measure for non-Gaussian fluctuations of the order parameter. At equilibrium in the thermodynamic limit it acquires two different universal values in the two phases: The Gaussian value 0 in the disordered phase and the value 2/3 in the ordered phase. At the transition point the Binder cumulant is scale invariant and it is a very convenient numerical probe for the existence of an equilibrium transition [28]. We will show that also in this non-equilibrium context for the 2d short-range case it behaves in the same way and allows to probe the existence of a transition.

III. DISCRETE TRUNCATED WIGNER APPROXIMATION
Before getting into the discussion of the results, it is useful to recap the basic ideas behind the DTWA and to discuss the accuracy of this method for this problem. In the following, we first review methodological details of the DTWA and afterwards we use exact diagonalization and MPS-TDVP data for a quantitative comparison.

A. DTWA Method
The DTWA is a semiclassical approximation, which has been used in many contexts concerning long-range interacting spin systems and has given noteworthy results, in terms of comparison with exact results and scalability to large system sizes. Precise details on the background can be found in [26,27], here we outline our concrete implementation. All the analysis is based on the construction of the discrete Wigner representation [25] which is a generalization to discrete Hilbert space of the usual Wigner representation (details can be found in Ref. 29). Summarizing, Wooters has shown that, given a discrete Hilbert space, the quantum dynamics can be represented through a discrete basis of operators. In the case of a single 1/2 spin, a possible basis choice iŝ where s β can take the values 1 1 1 , −1 1 −1 , 1 −1 −1 and −1 −1 1 andσ = σ xσyσz . With this basis choice, the expectation of any operator O acting on the Hilbert space of the single spin can be written as where w β ≡ 1 2 Tr Â βρ is the Wigner function, O w β (t) = Tr Â βÔ (t) are the Weyl symbols andÔ(t) ≡ e iĤtÔ e −iĤt . This representation can be extended also to our case of N spins considering as basis operatorŝ and writing as before the expectation of any operatorÔ acting on the Hilbert space of the N spins as Up to now everything is exact. The DTWA amounts to approximate the time-evolved basis operators as factorized objectŝ (13) The s ν βj (t) are initialized with the value for the corresponding β j given in Eq. (8) and obey a simple classical Hamiltonian dynamics given bẏ Here the symbol {· · · , · · · } is the Poisson bracket, µνρ is the Levi-Civita fully antisymmetric tensor, the variables s µ j, βj obey the angular-momentum Poisson brackets {s µ j, βj , s ν l, β l } = δ j l µνρ s ρ j, βj and the classical effective Hamiltonian is defined as For instance, the total longitudinal magnetization Eq. (5) can be evaluated in the DTWA scheme as In this form it still unpractical from the numerical point of view because the index β runs over 4 N values, so the sum would be unfeasible for large system sizes. The solution comes from the fact that the relation β w β = 1, so, in the cases when w β ≥ 0, it behaves as a probability distribution and it can be sampled through Monte Carlo sampling. With the initialization we choose we are in one of these lucky cases (see [26] for more details) and we can write Eq. (16) as the average over n r random initializations where each s j, βj is initialized with probability 1/2 in the condition 1 1 1 and probability 1/2 in the condition 1 −1 −1 .
Remarkably, the error bars do not scale with the system size, so this method is feasible also in the case of large systems. Moreover, results converge with a small number of randomness realizations (n r ); we show an example of this convergence in Appendix A. Unless otherwise specified here we use n r = 504.
Finally different sampling schemes, related to different choices of the operator in Eq. (8), can be employed. In Appendix A we briefly discuss these possible choices. All the results presented in the paper are essentially independent on the sampling method. Unless we specify otherwise, throughout the paper we use the sampling scheme specified in Eq. (8).
In the following we are going to compare the DTWA method with the results of other numerical methods in order to show its value also in our case.

B. Comparison with other methods
The comparison was done only in the case of onedimensional power-law interaction. In this case, in addition to the possibility to have results from exact diagonalization, it is possible to compare our data with tensornetwork (the MPS-TDVP) results [10,30] for larger sizes.
First of all we consider the case α = 0 of infiniterange interactions. In this limit the model reduces to the Lipkin-Meshkov-Glick model whose exact diagonalization dynamics can be easily studied. Being all the siteexchange operators conserved, there is an over-extensive number of constants of motion and the dynamics is integrable. Thanks to the conservation of the modulus of the total spin, the quantum dynamics is restricted to a Hilbert subspace whose dimension scales linearly with the system size and then large system sizes are feasible (for a detailed explanation see for instance [26]). We show some instances of comparison in Fig. 1. Let's first consider the case N = 100. We see that the curves of m x (t) deviate quite soon from each other, both for h < 1 and h > 1, but the time average (the one we are interested in) is actually the same.
We show also results for N = 10. Here we can see that in the ED case a phenomenon appears which is not captured by DTWA, the Rabi oscillations. Indeed, in this system an extensive number of eigenstates breaks the Z 2 symmetry in the thermodynamic limit. For any finite size, the true eigenstates are the even and odd superposition of these symmetry-breaking states and are separated by an exponentially small gap. Preparing the system in a symmetry-breaking state (as the one in Eq. (4)) gives rise therefore to Rabi oscillations of the magnetization with a frequency equal to the gap. Being this gap exponentially small in the system size, we cannot see these oscillations in panel (a), where the size is N = 100 and the gap is fantastically small. But we can see them in panel (c) and they are not caught by DTWA.
The existence of the Rabi oscillations is intimately related to the existence of a Z 2 symmetry, and the presence of resonant symmetry-breaking states put in interaction by the term with the h field. Explicitly breaking the symmetry breaks the resonance and there are no more oscillations. We do this in panel (c) where we show also a curve of m x (t) obtained adding to the Hamiltonian a small symmetry breaking term δh jσ x j . We see that there are no Rabi oscillations and the comparison with DTWA in terms of average is very good. So, in some sense, in DTWA one implicitly adds to the Hamiltonian a small symmetry breaking term. This is just what we operatively do when we want to see a quantum phase transition. We add a small symmetry-breaking term, we go to the thermodynamic limit, and then we send the small symmetry-breaking term to 0. Because we are interested here in the existence of a dynamical quantum phase transition with Z 2 symmetry breaking, this is exactly what we should do. DTWA does this implicitly for us, and in the thermodynamic limit the presence of a small symmetry-breaking term makes no difference both for DTWA and ED.
For α = 0 we can compare our DTWA results for the transverse magnetization m x (t) with the corresponding ones obtained through the TDVP method [10,30] (see Fig.2) for the case of N = 100 sites. The time scales we consider are much shorter than the times exponential in N needed for seeing the Rabi oscillations. Let us start focusing on the case α = 1.5 [ Fig. 2(a)]. We see that in this case the agreement is quite good both inside the symmetry-breaking phase (h = 0.7) and outside it (h = 1.5). On the opposite, for α = 3 [ Fig. 2(b)] the agreement is very good only when h = 1.5. When h = 0.7, the DTWA result decays much more slowly than TDVP. The two methods are in agreement for small values of α, as we expected from the existing literature on DTWA. In order to show the very good agreement when α is small, we plot in Fig. 3 the time-averaged longitudinal magnetization m x versus h for α = 0.1 and α = 1.5 obtained through the two methods. In both cases we see a very good agreement between the two methods. So, in the small-α regime we are interested in, the DTWA compares very well with the known results obtained through TDVP. This gives us an opportunity, because while TDVP can be used for at most N = 200 (see [10]), DTWA can be pushed up to much larger sizes, thus offering an opportunity of an accurate finite-size scaling.
We conclude this section by comparing DTWA results for the two-dimensional short-range case with the dynamics obtained by means of artificial neural networks (ANN) [31]. We show an example of comparison in Fig. 4 with data taken from Ref. [31] at large transverse fields in the regime where the Ising symmetry is restored in the long-time limit. As one can see, the DTWA compares well with the numerically exact ANN data.

IV. RESULTS
In this Section we will illustrate our results for the DPT obtained through the DTWA. We first analyze the one-dimensional long-range case [see Eq. (2)]. Later we will analyze the two-dimensional case with short-range interaction, Eq. (3). In this second case, we use also the Binder cumulant to get more reliable indications of the DPT. In both cases we address the steady state properties, and consider the behaviour of the time-averaged magnetization (6). We consider averages over a time T such that the magnetization has already converged and we specify it in any of the considered cases, explicitly studying the convergence in T for α 2.

A. Long-range model in 1D
Let us first analyze the one-dimensional long-range case and study the finite-size scaling of m x as a function of the transverse field h. We start with the case α = 0, where we can compare DTWA with the exact solution.
In form of the type The possible value β ∼ 0 implies logarithmic corrections of the form m x, N (h) ∼ (1/ log N )f (·). We obtain the scaling collapse shown in Fig. 6. In accordance with the exact solution, the DTWA reproduces the logarithmic corrections (see Fig. 6). Furthermore, we get a scaling exponent δ = 0.47 in good agreement with the exact exponent δ = 0.5 and a critical field h c = 1 corresponding to the exact result. The data collapse, shown in Fig. 6(b), is excellent. For comparison the same scaling is shown for the exact solution in Fig. 6(c). In the following of this section we are going to apply the methods illustrated here to the case with α = 0.
We first focus on the values of α ≤ 1. We show some examples of m x versus h for different sizes N and different α in Fig. 7 (also in this case the data shown are obtained for T = 200 where the observables have already attained their stationary value). Before doing the finite-size scal-ing, let us discuss more qualitatively what happens. For α = 0.1 [ Fig. 7 (a,c)] and α = 0.5 [ Fig. 7(b,d)] we observe a behaviour very similar to the case α = 0 shown in Fig. 5. In both cases the curves show a crossing at h c ∼ 1, the mean-field value. The tiny deviations from the mean-field are not relevant, only due to the fitting procedure. Indeed we can perform a finite-size scaling with the same method used for α = 0 and with the same scaling function as in Eq. (17). In the same Fig. 7 (lower panel) we show the collapsed curves. For α 0.5 the critical behaviour is mean-field like.
A different behaviour is observed at larger α (shorterrange interactions). We show the data for α = 1 in Fig. 8. The crossing point is clearly visible albeit the quality of the data collapse is not as good as in the previous cases. Several points are worth to be discussed. First of all the crossing field is still very close to one. The exponent δ ∼ 0.76, however, deviates significantly from the mean-field value. Although the DTWA does not allow ascertain how sizeable is the error in the scaling analysis, one can be confident in stating that for these parameters α there is still a transition point but the critical behaviour deviates from the mean-field.
Another feature that is worth noticing is that there is a range of transverse fields (in the disordered region above the critical field) where the magnetization becomes negative. Our analysis cannot exclude that this "reentrant" behaviour might still be a feature of the DTWA approximation, not present in more accurate analysis. It is however to be noted that this overshooting of the magnetization might be reminiscent of the chaotic behaviour observed in the mean-field dynamics of this model [32].
We conclude the analysis of the one-dimensional model by discussing the case of shorter-range interactions, α 2. In this regime our DTWA analysis indicates that the transition disappears. What we observe is that the (hypothetical) value of h at which the ordered phase disappears drops rapidly to zero. From the numerics, obviously, we cannot exclude that a transition point still exists at h c 1. As we will discuss in the following, this scenario is however quite unlikely.
In order to show that in the steady state the magnetization vanishes for (ideally) any value of h, we show in Fig. 9 the dependence of the magnetization as a function of the averaging time T of Eq.(6), for several different values of h. We consider two prototypical cases, α = 2 [panel (a)] and α = 3 [panel (b]. The data clearly show that in the limit T → ∞ the magnetization drops to zero. The slope decreases with h and at small value of the field one should go to larger times T to see this trend. In the insets a zoom of the curves at h = 0.3 [panel (a)] and h = 0. 2 [panel (b] clearly indicates the trend. The data, shown for a system of N = 100, are quite generic in this regime (finite size corrections are small). In order to push the bound of the ordered phase to smaller field one should consider larger systems and larger T . As previously mentioned, in principle we cannot exclude a critical point at h ∼ 0.1. We find this case, however, improbable, see however Ref. [33]. To further confirm this picture we analyzed also the case of a one-dimensional system with nearest-neighbour interaction where no transition exists. The results are strictly similar to those of Fig. 9 and are summarized in Appendix B. While in 1d the short-range case is trivial (no DPT is expected), the picture changes drastically by moving to higher dimensions. In the next Section we consider the case of a two-dimensional short-range interacting system as defined in Eq. (3).

B. Two-dimensional short-range model
This case is of particular importance for several reasons. First of all, to our knowledge, it has never be considered so far. Moreover, we expect that the transition will deviate from the mean-field behaviour. This then leads to the question whether the DTWA is capable to detect the transition and its non-mean field type character. If this is the case, a very important question to understand is if the system thermalizes and the dynamical transition corresponds to a thermal-equilibrium transition. The discussion below will try to address some of these points by analysing both the magnetization and the Binder cumulant.
In Fig. 10 we show the behaviour of the magnetization as a function of 1/L for different values of transverse field. Here we take T = 6 · 10 5 due to the long convergence times (this is essentially the limiting factor that forbids us to consider larger lattice sizes). DTWA indicates the existence of a transition for h * 0.7. In the ordered phase the magnetization increases with the system size and tends to converge only for the largest samples. This type of finite-size effects were observed also in the onedimensional case where the convergence with size was similarly attained only for N ∼ 100 − 200.
The value of the critical field is not far from the value obtained with quantum Monte Carlo simulations at thermal equilibrium. In this framework, taking a temperature T (h) such that the thermal energy coincides with the value of the energy in the DTWA dynamics, we see a transition between ordered and disordered phase at h *

38].) The equilibrium Binder cumulant is defined as
where · · · T (h) is the thermal-equilibrium average at the temperature T (h) defined above. The curves at different system size cross each other around h = h * Th . The agreement between dynamical and thermal transition suggests that the model is quantum chaotic and thermalizing. We can show the presence of quantum chaos by considering the level spacing distribution and checking that it is near to the Wigner-Dyson one [39]. We can do that by means of a probe called average level spacing ratio r (see [40] for definition and discussion). Using exact diagonalization in the fully symmetric Hilbert subspace of a 5 × 4 model we find r 0.52, a value very near to the quantum-chaotic Wigner-Dyson one r WD = 0.5295. It is therefore not surprising that the quantum dynamics shows a transition closely corresponding to the thermal one.
Together with the magnetization, we also examine the time-averaged Binder cumulant U L , defined in Eq. (7), and we plot it versus h. We have already discussed the thermal-equilibrium quantum Monte Carlo results in Fig. 11. We show data for DTWA in Fig. 12 and the ones for exact diagonalization in Fig. 14, considering in both cases different system sizes. Let us first focus on Fig. 12. The curves cross around h * 0.6 and the Binder cumulant rapidly decreases with L at larger fields. This is physically sound: The total magnetization is the sum of the local magnetizations which behave as uncorrelated random variables at large h because the correlation length is very short. The sum of uncorrelated random variables tends to a Gaussian as the number of random The Binder cumulant has been evaluated averaging over a time (T = 10 4 ) shorter than the time needed to attain an asymptotic value in the DTWA scheme. The point is that, before this asymptotic value, the DTWA Binder cumulant attains a metastable plateau: We show some examples in Fig. 13. This plateau gives rise to the crossing behaviour we can see in Fig. 12 while the asymptotic value does not. The metastable plateau therefore shows a behaviour more similar to the ones given by quantum Monte Carlo (Fig. 11) and by exact diagonalization (Fig. 14). This suggests that in this context DTWA gives physically more sound results for a finite time, although the difference between the metastable plateau and the asymptotic value is very small.
We study the behaviour of the Binder cumulant also by means of exact diagonalization. In Fig. 14 we show the exact-diagonalization Binder cumulant versus h for small system sizes. The trend is the same as that observed in Fig. 12. The point where there is a crossing between curves for different system sizes occurs at h * ∼ 0.6, which is in good agreement with the value found using DTWA. We stress again that for increasing system size the Binder cumulant tends to 2/3 in the ordered phase and to 0 in the disordered one, exactly as it occurs in the thermalequilibrium case.
In conclusion, for the two-dimensional short-range case there is a dynamical transition closely corresponding to the thermal one due to the fact that the system appears quantum chaotic and thermalizing. Remarkably, DTWA can see the existence of this transition.

V. CONCLUSIONS
In conclusion we have used DTWA to study the dynamical quantum phase transition in Ising spin models. Our aim was exploring the existence of a transition between an ordered and a disordered phase in the steady state and the properties of this transition focusing on a local order parameter, the time-averaged longitudinal magnetization.
We have first focused on the long-range onedimensional case where interactions decay with the power α of the distance. Here we have compared DTWA with numerically exact results (exact diagonalization for α = 0 and TDVP) and we have found a good agreement. Thanks to the good scalability of DTWA, we have done a finite-size scaling of the time-averaged longitudinal magnetization and we have studied the critical exponents of the transition between ordered and disordered phase. For α small (α = 0.1, 0.5) we have found the same critical exponents as the mean-field case (α = 0). For α = 1 we have found critical exponents significantly different from the mean-field case and we have found that the magnetization changes sign in the critical region. We do not know if this is a physical result or an effect of the DTWA approximation which should not work very well in the critical region due to the long correlation length of the physical system. For α 2 we have found no scaling at all with the system size and we have put a lower bound to the value of h for which the longitudinal magnetization vanishes at long times. We argue that this is most probably the case also for smaller h but we cannot see it due to the extremely long convergence times in the DTWA scheme (this is the same situation occurring if we apply DTWA to a short range one-dimensional Ising model).
We further considered the 2d short-range model, not considered in this context so far, again applying the DTWA approximation. Our data confirms that the DTWA is able to capture the existence of a transition and the value of the critical field compares well with the one of a corresponding thermal transition. We argued that this is physically sound showing that the model is quantum chaotic by means of exact diagonalization. In order to attempt a scaling analysis and thus to confirm that the associated critical exponents are the thermal ones it would be necessary to consider even larger systems, which might be an interesting prospect for the future. As discussed in Sec. III, in the DTWA approach one has to solve classical equations of motions for different random initial configuration. Physical quantities are obtained upon averaging over this initial distribution. In this Appendix we report on some details of the sampling procedure we used to obtain the results reported in the body of the paper.
First of all it is important to understand how the results depend on the number of random initial realizations n r . In Fig. 15 we consider the dependence of the average magnetization as a function of the number of initializations n r . We show the case of α = 0.1; the behaviour is however quite generic. Away from the critical field h c , the order parameter m x converges very rapidly to its asymptotic value and no significant changes happen by increasing n r . Since we are interested in determining transition points, the behavior m x as a function of n r is more notable in the critical region. Close to the transition point the convergence with the number of realizations is slower. In any case after few hundreds initial configurations the results seem stable. We choose n r = 504 for most of the calculations, if not stated otherwise.
In addition to the number of initial configurations over which performing the sampling, another aspect to con- sider is the choice of the sampling scheme. Indeed, using phase point operatorÂ α one can map each basis state of Hilbert space to a point in phase space. There are different possible choice of this phase operator and the one shown in Eq. (8) is not the only one. Any other possible choice for phase operator can be derived by some unitary transformation,Â β =ÛÂ βÛ † . In [41] the following phase operator was considered (more details about this construction can be found there) where s β can take the values 1 −1 1 , −1 −1 −1 , 1 1 −1 and −1 1 1 and σ = σ xσyσz which is obtained by flipping the sign of the second component of s β . Fig.16 show the comparison, as a function of the averaging time T for two different values of α. DTWA is further compared to exact diagonalization. In the fully connected case (α = 0) the different samplings lead to essentially the same result and agree with the diagonalization data. Smaller distances are observed in the bottom panel for the case α = 1. It should be noted that deviations appear only at the the third decimal digit. These differences may be important only very close to the transition point and may also contribute to the uncertainties in the scaling plots that we observe for α ∼ 1. However, In the case of spin chain with short-range interaction we do not expect to have ordered non-equilibrium steadystate (it corresponds to the one-dimensional model studied in the paper when α is very large). It is however useful to check this result with DTWA as an additional test of its quality. Following the same approach used to argue the absence of a critical point for α 2 we analyse how the magnetization scales with T for different values of the transverse field. The result of this analysis is presented in Fig. 17. Down to h = 0.1 the steady state magnetization (at large T ) tends to zero (top panel). The inset in the top panel shows that in order to see the suppression of the magnetization at large T one should go to very large values. In the top panel we considered a chain of length N = 100. Because of the short-range correlations in this case, the behaviour is essentially independent on N as displayed by the bottom panel of Fig. 17.