Disentanglement approach to quantum spin ground states: field theory and stochastic simulation

Consider a quantum spin system with Hamiltonian $\hat{H}$. The ground state |ψ_G⟩ of the system can be obtained from a generic state |ψ₀⟩ by performing imaginary time evolution: since at late imaginary times τ all excited states are exponentially suppressed compared to the ground state, one has

$$\begin{equation}\vert {\psi }_{\mathrm{G}}\rangle \sim \underset{\tau \to \infty }{\mathrm{lim}}\enspace {\mathrm{e}}^{-\hat{H}\tau }\vert {\psi }_{0}\rangle ,\end{equation} \tag{ 1 }$$

where we set ℏ = 1. It is then natural to introduce the Euclidean time evolution operator $\hat{U}\left(\tau \right)={\mathrm{e}}^{-\hat{H}\tau }$. Without loss of generality, let us consider initial states |ψ₀⟩ that are product states; these can be conveniently parameterized in terms of a single reference state |⇓⟩ ≡ ∏_i |↓⟩_i , where ${\hat{S}}_{i}^{-}\vert {{\downarrow}\rangle }_{i}=0$. The Euclidean time evolution from an arbitrary state |ψ₀⟩ is then obtained by considering the modified time evolution operator

$$\begin{equation}\hat{\mathcal{U}}\left(\tau \right)\equiv \hat{U}\left(\tau \right){\hat{U}}_{0},\end{equation} \tag{ 2 }$$

where the unitary operator ${\hat{U}}_{0}$ satisfies

$$\begin{equation}{\hat{U}}_{0}\vert {\Downarrow}\rangle =\vert {\psi }_{0}\rangle .\end{equation} \tag{ 3 }$$

The ground state expectation value of an observable $\hat{\mathcal{O}}$ can thus be written as

$$\begin{equation}{\mathcal{O}}_{\mathrm{G}}=\underset{\tau \to \infty }{\mathrm{lim}}\enspace \frac{\langle {\Downarrow}\vert {\hat{\mathcal{U}}}^{{\dagger}}\left(\tau \right)\hat{\mathcal{O}}\hat{\mathcal{U}}\left(\tau \right)\vert {\Downarrow}\rangle }{\langle {\Downarrow}\vert {\hat{\mathcal{U}}}^{{\dagger}}\left(\tau \right)\hat{\mathcal{U}}\left(\tau \right)\vert {\Downarrow}\rangle }.\end{equation} \tag{ 4 }$$

The denominator of equation (4) provides the necessary normalization, since $\hat{\mathcal{U}}\left(\tau \right)$ inherits the non-unitarity of $\hat{U}\left(\tau \right)$. All information about ground state expectation values is then encoded in the late-time behavior of $\hat{\mathcal{U}}\left(\tau \right)$. Let us consider spin systems with a quadratic Hamiltonian

$$\begin{equation}\hat{H}=-J\sum _{ijab}{\mathcal{J}}_{ij}^{ab}{\hat{S}}_{i}^{a}{\hat{S}}_{j}^{b}-\sum _{ia}{h}_{i}^{a}{\hat{S}}_{i}^{a},\end{equation} \tag{ 5 }$$

where a, b ∈ {+, −, z} and the spin operators ${\hat{S}}_{j}^{a}$ on site j satisfy the SU(2) commutation relations $\left[{\hat{S}}_{j}^{z},{\hat{S}}_{{j}^{\prime }}^{{\pm}}\right]={\pm}{\delta }_{j{j}^{\prime }}{\hat{S}}_{j}^{{\pm}}$, $\left[{\hat{S}}_{j}^{+},{\hat{S}}_{{j}^{\prime }}^{-}\right]=2{\delta }_{j{j}^{\prime }}{\hat{S}}_{j}^{z}$. We consider a symmetric interaction matrix ${\mathcal{J}}_{ij}^{ab}$ with interaction strength J and an applied magnetic field ${h}_{i}^{a}$. For systems whose Hamiltonian is of the form (5), the operator $\hat{\mathcal{U}}\left(\tau \right)$ can be conveniently re-expressed using a disentanglement formalism [24, 55, 56], recently applied in the context of real time evolution [58, 59]. By performing a Hubbard–Stratonovich (HS) decoupling [60, 61] followed by a Lie-algebraic disentanglement transformation [24, 55, 56, 62, 63], $\hat{\mathcal{U}}\left(\tau \right)$ can be exactly represented as a functional integral [24, 55, 56]:

$$\begin{equation}\hat{\mathcal{U}}\left(\tau \right)=\int \mathcal{D}\varphi \enspace {\mathrm{e}}^{-{S}_{0}\left[\varphi \right]}{\hat{\mathcal{U}}}^{\mathrm{s}}\left(\tau \right),\end{equation} \tag{ 6 }$$

where the noise action S₀ is given by ¹

$$\begin{equation}{S}_{0}\left[\varphi \right]\equiv \frac{J}{4}{\int }_{0}^{\tau }\mathrm{d}{\tau }^{\prime }\sum _{abij}{\left({\mathcal{J}}^{-1}\right)}_{ij}^{ab}{\varphi }_{i}^{a}\left({\tau }^{\prime }\right){\varphi }_{j}^{b}\left({\tau }^{\prime }\right)\end{equation} \tag{ 7 }$$

and the stochastic time evolution operator ${\hat{\mathcal{U}}}^{\mathrm{s}}$ is defined as a product of on-site operators:

$$\begin{equation}{\hat{\mathcal{U}}}^{\mathrm{s}}\left(\tau \right)\equiv \prod _{j}{\hat{\mathcal{U}}}_{j}^{\mathrm{s}}\left(\tau \right)=\prod _{j}{\text{e}}^{{\xi }_{j}^{+}\left(\tau \right){\hat{S}}_{j}^{+}}\enspace {\text{e}}^{{\xi }_{j}^{z}\left(\tau \right){\hat{S}}_{j}^{z}}\enspace {\text{e}}^{{\xi }_{j}^{-}\left(\tau \right){\hat{S}}_{j}^{-}}.\end{equation} \tag{ 8 }$$

The operators ${\hat{\mathcal{U}}}_{j}^{\mathrm{s}}$ have a functional dependence on the fields $\varphi =\left\{{\varphi }_{i}^{a}\right\}$ via the disentangling variables $\xi \equiv \left\{{\xi }_{j}^{\nu }\right\}$, which satisfy [24]

$$\begin{equation}{\dot {\xi }}_{j}^{+}={{\Phi}}_{j}^{+}+{{\Phi}}_{j}^{z}{\xi }_{j}^{+}-{{\Phi}}_{i}^{-}{{\xi }_{j}^{+}}^{2},\end{equation} \tag{ 9a }$$

$$\begin{equation}{\dot {\xi }}_{j}^{z}={{\Phi}}_{j}^{z}-2{{\Phi}}_{j}^{-}{\xi }_{j}^{+},\end{equation} \tag{ 9b }$$

$$\begin{equation}{\dot {\xi }}_{j}^{-}={{\Phi}}_{j}^{-}\enspace \mathrm{exp}\enspace {\xi }_{j}^{z},\end{equation} \tag{ 9c }$$

where ${{\Phi}}_{j}^{a}={h}_{j}^{a}+J{\varphi }_{j}^{a}$. The initial conditions of the disentangling variables are determined from $\hat{\mathcal{U}}\left(0\right)={\hat{U}}_{0}$; for example, for a spin-1/2 system, the general product state |ψ₀⟩ ≡ ∏_i (a_i , b_i ) corresponds to the initial conditions

$$\begin{equation}{\xi }_{i}^{+}\left(0\right)={a}_{i}/{b}_{i},\end{equation} \tag{ 10a }$$

$$\begin{equation}{\xi }_{i}^{z}\left(0\right)=-2\enspace \mathrm{log}\left({b}_{i}\right),\end{equation} \tag{ 10b }$$

$$\begin{equation}{\xi }^{-}\left(0\right)=-{a}_{i}^{{\ast}}/{b}_{i}^{{\ast}}.\end{equation} \tag{ 10c }$$

For completeness, we outline the derivation of equations (6) and (9) in appendix A. Equation (6) can be seen as an exact path integral representation of the time-evolution operator. The operators inside the functional average (6) are decoupled over sites and act in a simple way on any state of interest. While each trajectory parameterized by ${\xi }_{i}^{a}$ describes a non-interacting spin and has no entanglement, the effect of interactions is encoded in the correlations between the fields ${\varphi }_{i}^{a}$, determined by the noise action (7); the interacting quantum spin system is thus fully retrieved upon performing the functional integral in (6). Equation (6) allows one to formulate an exact field theoretical description of lattice spin systems, as we show in section 3. The noise action (7) can be diagonalized in terms of a new set of fields $\phi =\left\{{\phi }_{i}^{a}\right\}$ by performing the linear transformation ${\phi }_{i}^{a}={\sum }_{bj}{O}_{ij}^{ab}{\varphi }_{j}^{b}$, where O is a matrix satisfying ${O}^{\mathrm{T}}{\mathcal{J}}^{-1}O/2J=\mathbb{1}$. With this transformation, the noise action takes the form [58]

$$\begin{equation}{S}_{0}\left[\phi \right]\equiv \frac{1}{2}{\int }_{0}^{\tau }\mathrm{d}{\tau }^{\prime }\sum _{ai}{\phi }_{i}^{a}\left({\tau }^{\prime }\right){\phi }_{i}^{a}\left({\tau }^{\prime }\right).\end{equation} \tag{ 11 }$$

Notably, due to the Gaussian nature of (11), the fields ϕ can also be interpreted as delta-correlated, unit-variance Gaussian white noise variables [24, 55]. Thus, the functional integral in equation (6) can be equivalently seen as an average over the stochastic processes ${\phi }_{i}^{a}\left(\tau \right)$ [24, 55, 58]:

$$\begin{equation}\hat{\mathcal{U}}\left(\tau \right)={\left\langle \prod _{j}{\text{e}}^{{\xi }_{j}^{+}\left(\tau \right){\hat{S}}_{j}^{+}}\enspace {\text{e}}^{{\xi }_{j}^{z}\left(\tau \right){\hat{S}}_{j}^{z}}\enspace {\text{e}}^{{\xi }_{j}^{-}\left(\tau \right){\hat{S}}_{j}^{-}}\right\rangle }_{\phi }.\end{equation} \tag{ 12 }$$

The notation ⟨...⟩_ϕ denotes averaging with respect to the noise action (11), so that the relevant probability law is the Wiener measure. The equations of motion equation (9) are then interpreted as SDEs for the variables ${\xi }_{i}^{a}$ [24, 55]. By representing each of the time evolution operators in equation (4) using the disentangling formula (12), one can express quantum ground state expectation values as classical averages. Introducing independent sets of forwards and backwards fields (labeled by the additional subscripts f, b), ${\phi }_{\mathrm{f}}\equiv \left\{{\phi }_{\mathrm{f},i}^{a}\right\}$ and ${\phi }_{\mathrm{b}}\equiv \left\{{\phi }_{\mathrm{b},i}^{a}\right\}$, and the associated disentangling variables ${\xi }_{\mathrm{f}}\equiv \left\{{\xi }_{\mathrm{f},i}^{a}\right\}$, ${\xi }_{\mathrm{b}}\equiv \left\{{\xi }_{\mathrm{b},i}^{a}\right\}$, one has

$$\begin{equation}{\mathcal{O}}_{\mathrm{G}}=\underset{\tau \to \infty }{\mathrm{lim}}\enspace \frac{\langle {F}_{\mathcal{O}}{\left(\tau \right)\rangle }_{{\phi }_{\mathrm{f}},{\phi }_{\mathrm{b}}}}{\langle {F}_{\mathbb{1}}{\left(\tau \right)\rangle }_{{\phi }_{\mathrm{f}},{\phi }_{\mathrm{b}}}},\end{equation} \tag{ 13 }$$

where the classical function ${F}_{\mathcal{O}}$ corresponding to the operator $\hat{\mathcal{O}}$ is defined by

$$\begin{equation}{F}_{\mathcal{O}}\equiv \langle {\Downarrow}\vert {\left[{\hat{\mathcal{U}}}^{\mathrm{s}}\left(\tau \right)\right]}^{{\dagger}}\hat{\mathcal{O}}{\hat{\mathcal{U}}}^{\mathrm{s}}\left(\tau \right)\vert {\Downarrow}\rangle ,\end{equation} \tag{ 14 }$$

and ${F}_{\mathbb{1}}$ is obtained from (14) when $\hat{\mathcal{O}}$ is replaced by the identity operator. It can be readily seen that the functions ${F}_{\mathcal{O}}$ take the same functional form as their real time counterparts, given in [58, 59]. Here, in contrast to said references, we express all initial states in terms of a single reference state and variable initial conditions ${\xi }_{i}^{a}\left(0\right)$: in this way, there is a one-to-one correspondence between observables and their classical counterparts, regardless of the initial state. The reference state |⇓⟩ was selected because it results in the simplest classical expressions [58, 59]. For example, for spin-1/2 systems one has ${\hat{S}}_{i}^{a}={\sigma }_{i}^{a}/2$, where ${\sigma }_{i}^{a}$ are the Pauli matrices. In this case, by acting with (12) on the reference state |⇓⟩, the classical function for the normalization is found to be

$$\begin{equation}{F}_{\mathbb{1}}\left(\tau \right)=\prod _{i}\left[1+{\xi }_{\mathrm{f},i}^{+}\left(\tau \right){\xi }_{\mathrm{b},i}^{+{\ast}}\left(\tau \right)\right]{\mathrm{e}}^{-\frac{1}{2}\left[{\xi }_{\mathrm{f},i}^{z}\left(\tau \right)+{\xi }_{\mathrm{b},i}^{z{\ast}}\left(\tau \right)\right]},\end{equation} \tag{ 15 }$$

while the longitudinal magnetization ${\mathcal{M}}_{z}={\sum }_{j}\langle {\hat{S}}_{j}^{z}\rangle /N$ corresponds to the classical function [58]

$$\begin{equation}{F}_{{\mathcal{M}}_{z}}\left(\tau \right)=\frac{{F}_{\mathbb{1}}\left(\tau \right)}{N}\sum _{j}\frac{1-{\xi }_{\mathrm{f},j}^{+}\left(\tau \right){\xi }_{\mathrm{b},j}^{+{\ast}}\left(\tau \right)}{1+{\xi }_{\mathrm{f},j}^{+}\left(\tau \right){\xi }_{\mathrm{b},j}^{+{\ast}}\left(\tau \right)}.\end{equation} \tag{ 16 }$$

For any observable, the appropriate classical function can be constructed using the building blocks provided in [59]. The expectation values of functions such as (15) and (16) can be evaluated numerically by averaging them over realizations of the stochastic processes ${\phi }_{i}^{a}$, as done in [58, 59] for real time evolution; the processes ${\phi }_{i}^{a}$ determine the time evolution of the variables ${\xi }_{i}^{a}$ via the SDEs (9). The key toward both analytical and numerical developments is identifying the trajectories ϕ which provide the largest contribution to each observable; we discuss this in the following section.

In the disentanglement formalism, ground state expectation values can be numerically computed by averaging classical functions over realizations of trajectories ${\phi }_{i}^{a}$ distributed according to the action (11). In this approach, which we refer to as direct sampling, one preferentially generates trajectories that are close to the non-interacting limit ${\phi }_{i}^{a}\left(\tau \right)=0$, as shown in appendix B. This was done for real-time evolution in [58, 59]. However, the trajectories that contribute most significantly to a given functional integral may be substantially different from the non-interacting trajectories, as also pointed out in [64]. Therefore, going beyond previous applications of the disentanglement formalism [24, 55, 56, 58, 59], we presently identify the trajectories yielding the dominant contributions. To this end, it is convenient to work with the action (7) written in terms of the fields ${\varphi }_{i}^{a}$. An observable $\mathcal{O}$ can be written as

$$\begin{equation}\mathcal{O}\equiv \langle \hat{\mathcal{O}}\rangle =\langle {F}_{\mathcal{O}}{\left[\varphi \right]\rangle }_{\varphi }\end{equation} \tag{ 17 }$$

where φ = {φ_α } denotes all of the HS fields: the collective index α runs over sites, Lie algebra generators, and sets of fields (e.g. forwards and backwards). Equation (17) can be equivalently written as

$$\begin{equation}\mathcal{O}\equiv \int \mathcal{D}\varphi \enspace {\mathrm{e}}^{-{S}_{\mathcal{O}}\left[\varphi \right]},\end{equation} \tag{ 18 }$$

which defines the effective action ${S}_{\mathcal{O}}\left[\varphi \right]\equiv {S}_{0}\left[\varphi \right]-\mathrm{log}\enspace {f}_{\mathcal{O}}\left[\varphi \right]$ for the observable $\mathcal{O}$. The leading contribution to the integral (18) is given by the configuration φ_SP(τ) which extremizes the effective action:

$$\begin{equation}{\varphi }_{\text{SP}}:{\left.\frac{\delta {S}_{\mathcal{O}}}{\delta \varphi }\right\vert }_{{\varphi }_{\text{SP}}}=0.\end{equation} \tag{ 19 }$$

We refer to φ_SP as the SP field. The trajectory corresponding to the field φ_SP provides the leading order (LO) approximation to a given observable, as we further discuss in section 3.1.

One can regard equation (18) as an exact field theoretical formulation of the lattice quantum spin system; equation (18) can then be expanded about the SP to analytically obtain corrections to observables beyond LO, computed by using the associated Feynman rules and diagrams. In practice, the expansion can be carried out up to a desired order, providing approximate analytical results for ground state expectation values. This approach is discussed and explicitly applied in section 3.

Alternatively, one can view the disentanglement approach as a numerical tool. φ_SP can then be used to perform an exact change of variables in the functional integral (18): this amounts to preferentially sampling trajectories near φ_SP, which yield the largest contributions to the integral. Notably, in this approach one does not truncate to a given order in the fluctuations: the resulting equations are still exact and fully capture the corresponding quantum problem. The measure transformation amounts to an importance sampling scheme to improve the numerical efficiency of the stochastic approach, which is formally exact and whose practical accuracy is determined by the finite number of simulations one performs. The importance sampling method is discussed in section 4, where we show how this can be used to numerically access much larger system sizes than it is possible by sampling according to the naive measure (11).

In the next section, we show how the SP trajectory is computed; this provides the starting point for both analytical and numerical developments.

We wish to identify the leading contribution to a given functional integral by solving equation (19). For definiteness, we illustrate this by considering the quantum Ising model in D spatial dimensions, but the same formalism can be applied to any spin Hamiltonian of the form (5). The quantum Ising model is given by the Hamiltonian

$$\begin{equation}{\hat{H}}_{\text{I}}=-J\sum _{\langle ij\rangle }^{N}{\hat{S}}_{i}^{z}{\hat{S}}_{j}^{z}-{\Gamma}\sum _{j=1}^{N}{\hat{S}}_{j}^{x},\end{equation} \tag{ 20 }$$

where ⟨ij⟩ denotes nearest-neighbor interactions. We consider a system of N = N₁ ×⋯ × N_D spin-1/2 degrees of freedom on a D-dimensional hypercubic lattice, with periodic boundary conditions and ferromagnetic interactions J > 0. We begin by considering the integrable one-dimensional case, and then generalize our results to D > 1, where the model cannot be solved exactly. For D = 1, the model (20) reduces to the quantum Ising chain and is solvable in terms of free fermions [65]; this allows the exact computation of physical observables in the thermodynamic limit. In the present units, the Ising chain has a QPT at Γ = J/2. For this model, the general result (9) specializes to the Euclidean Ising SDEs [24, 58, 59]

$$\begin{equation}{\dot {\xi }}_{i}^{+}\left(\tau \right)=\frac{{\Gamma}}{2}\left(1-{{\xi }_{i}^{+}}^{2}\right)+J{\xi }_{i}^{+}{\varphi }_{i},\end{equation} \tag{ 21a }$$

$$\begin{equation}{\dot {\xi }}_{i}^{z}\left(\tau \right)=-{\Gamma}{\xi }_{i}^{+}+J{\varphi }_{i},\end{equation} \tag{ 21b }$$

$$\begin{equation}{\dot {\xi }}_{i}^{-}\left(\tau \right)=\frac{{\Gamma}}{2}\enspace \mathrm{exp}\enspace {\xi }_{i}^{z}.\end{equation} \tag{ 21c }$$

A natural choice of observable is the ground state energy density _G. This can be computed using equation (13), according to the general formalism outlined in section 2.1. Alternatively, _G can also be obtained as

$$\begin{equation}{{\epsilon}}_{\mathrm{G}}=-\underset{{\tau }_{f}\to \infty }{\mathrm{lim}}\enspace \frac{1}{N{\tau }_{f}}\enspace \mathrm{log}\enspace A\left({\tau }_{f}\right),\end{equation} \tag{ 22 }$$

where we defined the Euclidean Loschmidt amplitude A(τ_f ) for the initial state |ψ₀⟩ as

$$\begin{equation}A\left({\tau }_{f}\right)=\langle {\psi }_{0}\vert \hat{U}\left({\tau }_{f}\right)\vert {\psi }_{0}\rangle .\end{equation} \tag{ 23 }$$

By computing _G by means of equation (22), one only needs to consider a single time evolution operator, which corresponds to a single set of HS fields ${\phi }_{i}^{a}$. Thus, using equation (22) allows us to simplify the subsequent analytical developments for the purpose of illustrating the method. The same results can however be obtained from the general formalism of section 2.1, involving two time evolution operators; see appendix C. To further simplify our calculations, we choose the all-down initial state |ψ₀⟩ = ⊗_j |↓⟩_j ≡ |⇓⟩; in this case, the Loschmidt amplitude is given by the functional integral [58]

$$\begin{equation}A\left({\tau }_{f}\right)=\int \mathcal{D}\varphi \enspace {\mathrm{e}}^{-{S}_{0}\left[\varphi \right]}\enspace {\mathrm{e}}^{-\frac{1}{2}{\sum }_{i}{\xi }_{i}^{z}\left({\tau }_{f}\right)},\end{equation} \tag{ 24 }$$

where the equation of motion of ${\xi }_{j}^{z}$ is given by (21b) and the initial conditions are ${\xi }_{i}^{a}\left(0\right)=0$. Following the discussion of section 2.2, we write the Loschmidt amplitude (24) as

$$\begin{equation}A\left({\tau }_{f}\right)=\int \mathcal{D}\varphi \enspace {\mathrm{e}}^{-S\left[\varphi \right]}\end{equation} \tag{ 25 }$$

which defines the Euclidean Loschmidt action:

$$\begin{equation}S\left[\varphi \right]=\frac{J}{2}{\int }_{0}^{{\tau }_{f}}\mathrm{d}\tau \left[\frac{1}{2}\sum _{ij}{\left({\mathcal{J}}^{-1}\right)}_{ij}{\varphi }_{i}\left(\tau \right){\varphi }_{j}\left(\tau \right)-\frac{{\Gamma}}{J}\sum _{i}{\xi }_{i}^{+}\left(\tau \right)+\sum _{i}{\varphi }_{i}\left(\tau \right)\right].\end{equation} \tag{ 26 }$$

The variables ${\xi }_{i}^{+}$ featured in the action (26) are themselves functionals of φ_i , as determined by (21). It follows that S[φ] cannot be written in terms of a Lagrangian involving only the fields φ and their time-derivatives, and it is thus not possible to obtain Euler–Lagrange equations in the standard way. Rather, in order to obtain the SP field configuration, one directly extremizes the action (26) with respect to varying the field φ_i . This yields the Loschmidt SP equation

$$\begin{equation}{\varphi }_{i}\left({\tau }^{\prime }\right){\vert }_{\text{SP}}=\frac{{\Gamma}}{J}\sum _{j}{\mathcal{J}}_{ij}{\int }_{0}^{{\tau }_{f}}\mathrm{d}\tau {\left.\frac{\delta {\xi }_{j}^{+}\left(\tau \right)}{\delta {\varphi }_{j}\left({\tau }^{\prime }\right)}\right\vert }_{\text{SP}}-1,\end{equation} \tag{ 27 }$$

where we used ${\sum }_{j}\;{\mathcal{J}}_{ij}=1\enspace \enspace \forall \enspace i$ and $\delta {\xi }_{i}^{+}/\delta {\varphi }_{j}\propto {\delta }_{ij}$. The subscript SP denotes quantities that are evaluated at the SP. By varying equation (21a) with respect to φ_i (τ'), one obtains

$$\begin{equation}\frac{\delta {\xi }_{i}^{+}\left(\tau \right)}{\delta {\varphi }_{i}\left({\tau }^{\prime }\right)}=J{\xi }_{i}\left({\tau }^{\prime }\right)\theta \left(\tau -{\tau }^{\prime }\right){\mathrm{e}}^{-{\int }_{{\tau }^{\prime }}^{\tau }{\gamma }_{i}\left(s\right)\mathrm{d}s}\equiv {{\Xi}}_{i}\left(\tau ,{\tau }^{\prime }\right),\end{equation} \tag{ 28 }$$

where θ(τ) is the Heaviside step function and we defined

$$\begin{equation}{\gamma }_{i}\left(s\right)\equiv {\Gamma}{\xi }_{i}^{+}\left(s\right)-J{\varphi }_{i}\left(s\right).\end{equation} \tag{ 29 }$$

Due to the translational invariance of the model (20) and of the chosen initial state, at the SP all ${\xi }_{i}^{+}$ take the same value, ${\xi }_{i}^{+}{\vert }_{\text{SP}}\equiv {\xi }_{\text{SP}}^{+}$. From the translational symmetry of equation (27), it also follows that φ_i |_SP = φ_SP and Ξ_i |_SP = Ξ_SP. Hence, in the translationally invariant case the SP equation for the field φ_SP simplifies to

$$\begin{equation}{\varphi }_{\text{SP}}\left({\tau }^{\prime }\right)=\frac{{\Gamma}}{J}{\int }_{0}^{{\tau }_{f}}{{\Xi}}_{\text{SP}}\left(\tau ,{\tau }^{\prime }\right)\mathrm{d}\tau -1.\end{equation} \tag{ 30 }$$

From equations (30) and (28) one immediately obtains the boundary condition φ_SP(τ_f ) = −1; equations (30), (28) and (21) further imply that φ_SP(τ') must remain real-valued at all times. To the best of our knowledge, the functional equation (30) cannot be solved analytically. However, for the computation of ground states one is only interested in the limit τ_f → ∞. A recursive numerical solution of equation (30) shows that in this limit, away from a transient near τ = 0 and a boundary region at τ ≲ τ_f , the SP equation is dominated by a time-independent plateau value ϕ_P; see appendix C. We may then assume that late-time plateau values, denoted by a subscript P, dominate the integrals, and approximate

$$\begin{equation}{\int }_{0}^{{\tau }_{f}}{{\Xi}}_{\mathrm{P}}\left(\tau ,{\tau }^{\prime }\right)\approx {\int }_{0}^{{\tau }_{f}}\theta \left(\tau -{\tau }^{\prime }\right)J{\xi }_{\mathrm{P}}^{+}\enspace {\text{e}}^{-{\gamma }_{\mathrm{P}}\left(\tau -{\tau }^{\prime }\right)}=\frac{J{\xi }_{\mathrm{P}}^{+}}{{\gamma }_{\mathrm{P}}},\end{equation} \tag{ 31 }$$

where ${\gamma }_{\mathrm{P}}\equiv {\Gamma}{\xi }_{\mathrm{P}}^{+}-J{\varphi }_{\mathrm{P}}$. Convergence of the integral in equation (31) requires γ_P < 0. Assuming that this condition is satisfied, which can be self-consistently verified a posteriori, in the τ_f → ∞ limit equation (30) is reduced to an algebraic equation for φ_P:

$$\begin{equation}{\varphi }_{\mathrm{P}}=\frac{{\Gamma}{\xi }_{\mathrm{P}}}{{\Gamma}{\xi }_{\mathrm{P}}^{+}-J{\varphi }_{\mathrm{P}}}-1.\end{equation} \tag{ 32 }$$

In the absence of translational invariance, an analogous set of equations for the plateau fields can be obtained from (27). Equation (32) can be solved together with the condition that ${\xi }_{\mathrm{P}}^{+}$ is a fixed point of the Euclidean dynamics when φ = φ_P, yielding four solutions:

$$\begin{equation}\begin{aligned}\hfill & {\xi }_{\mathrm{P}}^{+}=\begin{cases}\left[J-\sqrt{{J}^{2}-{{\Gamma}}^{2}}\right]/{\Gamma}\quad \hfill \\ \left[J+\sqrt{{J}^{2}-{{\Gamma}}^{2}}\right]/{\Gamma}\quad \hfill \\ 1\quad \hfill \\ -1\quad \hfill \end{cases}\quad {\varphi }_{\mathrm{P}}=\begin{cases}-\sqrt{{J}^{2}-{{\Gamma}}^{2}}/J\quad \hfill \\ \sqrt{{J}^{2}-{{\Gamma}}^{2}}/J\quad \hfill \\ 0\quad \hfill \\ 0\quad \hfill \end{cases}\hfill \\ \hfill & {\gamma }_{\mathrm{P}}=\begin{cases}J\quad \hfill \\ J\quad \hfill \\ {\Gamma}\quad \hfill \\ -{\Gamma}.\quad \hfill \end{cases}\hfill \end{aligned}\end{equation} \tag{ 33 }$$

In order for the SP field to be real valued, the first and second solutions are only acceptable when Γ < J; they both give γ_P = J, and the corresponding ${\xi }_{\mathrm{P}}^{+}$ are reciprocal to each other. The fourth solution is not acceptable as it gives γ_P < 0: it corresponds to a maximum of the action (26). We refer to the first and second SPs as the small-Γ SPs and to the third one as the large-Γ SP, as they give the LO contribution to the ground state energy in these limits; see section 3.1 below. Notably, the plateau values φ_P in equation (33) coincide with the effective fields acting on each spin within the mean field (MF) approximation; similarly, the disentangling variables ${\xi }_{\mathrm{P}}^{+}$ parameterize the MF ground states. For comparison, we provide details of the MF solution in appendix D. This finding has a transparent physical interpretation: the path integral (24) is a sum over configurations of non-interacting spins, i.e. product states, and the SP trajectory is the single such configuration which gives the best approximation to the ground state energy. The product state which best approximates a ground state is precisely given by MF; the first and second SPs in (33), which have opposite φ_P, arise from spontaneous symmetry breaking at the MF level. This interpretation suggests that the correspondence between the SP configuration and MF is general. In order to compute ground state expectation values using the present method, it is therefore convenient to initialize the system in the MF ground state and subsequently perform imaginary time evolution toward the true ground state. This is tantamount to initializing the disentangling variables at their plateau values, ${\xi }_{i}^{a}\left(0\right)={\xi }_{\mathrm{P}}^{a}$, which removes the initial transient behavior. In principle, the above discussion should be repeated for every observable, since each corresponds to a different effective action and therefore to a different SP equation. However, it can be shown that the SP configuration of the Loschmidt action also extremizes the effective action for all physical observables, as obtained from the general formalism of section 2.1; see appendix C. The findings of this section also readily generalize to higher dimensions. The SP solution corresponds to MF also for D > 1; a detailed derivation is provided in appendix C. For instance, the plateau SP values for an isotropic quantum Ising model in D spatial dimensions are given by

$$\begin{equation}\begin{aligned}\hfill & {\xi }_{\mathrm{P}}^{+}=\begin{cases}\left[DJ-\sqrt{{D}^{2}{J}^{2}-{{\Gamma}}^{2}}\right]/{\Gamma}\quad \hfill \\ \left[DJ+\sqrt{{D}^{2}{J}^{2}-{{\Gamma}}^{2}}\right]/{\Gamma}\quad \hfill \\ 1\quad \hfill \\ -1\quad \hfill \end{cases}\hfill \\ \hfill & {\varphi }_{\mathrm{P}}=\begin{cases}-\sqrt{{D}^{2}{J}^{2}-{{\Gamma}}^{2}}/DJ\quad \hfill \\ \sqrt{{D}^{2}{J}^{2}-{{\Gamma}}^{2}}/DJ\quad \hfill \\ 0\quad \hfill \\ 0\quad \hfill \end{cases}\quad {\gamma }_{\mathrm{P}}=\begin{cases}DJ\quad \hfill \\ DJ\quad \hfill \\ {\Gamma}\quad \hfill \\ -{\Gamma}.\quad \hfill \end{cases}\hfill \end{aligned}\end{equation} \tag{ 34 }$$

We begin by considering the LO term; this is given by the plateau field configurations (33) obtained in section 2.3. For the remainder of this section, it is convenient to initialize the disentangling variables at their plateau values; since ${\xi }_{i}^{+}\left(0\right)={\xi }_{\mathrm{P}}^{+}$ corresponds to the MF ground state |MF⟩, this is equivalent to expressing the ground state energy density in the thermodynamic limit as

$$\begin{equation}{{\epsilon}}_{\mathrm{G}}=\underset{{\tau }_{f}\to \infty }{\mathrm{lim}}\enspace \underset{N\to \infty }{\mathrm{lim}}-\frac{1}{{\tau }_{f}N}\enspace \mathrm{log}\langle {\Downarrow}\vert \hat{U}\left({\tau }_{f}\right)\vert \text{MF}\rangle .\end{equation} \tag{ 35 }$$

Normalization of the initial state also implies ${\xi }_{i}^{z}\left(0\right)=\mathrm{log}\left(1+{\xi }_{\mathrm{P}}^{+2}\right)$. Equation (35) can thus be written as

$$\begin{equation}{{\epsilon}}_{\mathrm{G}}=\underset{{\tau }_{f}\to \infty }{\mathrm{lim}}\enspace \underset{N\to \infty }{\mathrm{lim}}-\frac{1}{{\tau }_{f}N}\enspace \mathrm{log}\enspace \mathcal{A}\left({\tau }_{f}\right)\end{equation} \tag{ 36 }$$

in terms of a modified Loschmidt amplitude, given by

$$\begin{equation}\mathcal{A}\left({\tau }_{f}\right)=\langle {\Downarrow}\vert \hat{\mathcal{U}}\left({\tau }_{f}\right)\vert {\Downarrow}\rangle .\end{equation} \tag{ 37 }$$

The modified time evolution operator $\hat{\mathcal{U}}\left({\tau }_{f}\right)$ in equation (37) is given by equation (2) with the condition $\hat{\mathcal{U}}\left(0\right)\vert {\Downarrow}\rangle =\vert \text{MF}\rangle$; due to this definition, the modified Loschmidt amplitude $\mathcal{A}\left({\tau }_{f}\right)$ corresponds to the same effective action (26) as A(τ_f ). The analysis of the action (26) in section 2.3 concerns the infinite time limit and is independent of the initial conditions. In this limit, we can again assume that all integrals are dominated by the plateau values; therefore, the earlier discussion equally applies to the present case, and the two amplitudes $\mathcal{A}\left({\tau }_{f}\right)$, A(τ_f ) are dominated by the same large-τ_f plateaus.

The LO approximation to the ground state energy density in the thermodynamic limit can then be obtained as

$$\begin{equation}{{\epsilon}}_{\mathrm{G}}\approx -\underset{{\tau }_{f}\to \infty }{\mathrm{lim}}\enspace \underset{N\to \infty }{\mathrm{lim}}\enspace \frac{1}{{\tau }_{f}N}\enspace \mathrm{log}\sum _{\text{s.p.}}{\mathrm{e}}^{-{S}_{\mathrm{P}}},\end{equation} \tag{ 38 }$$

where the sum runs over the different SPs and the plateau action is given by

$$\begin{equation}{S}_{\mathrm{P}}=\frac{NJ}{4}{\varphi }_{\mathrm{P}}^{2}{\tau }_{f}-\frac{{\Gamma}}{2}{\xi }_{\mathrm{P}}{\tau }_{f}+\frac{J}{2}{\varphi }_{\mathrm{P}}{\tau }_{f}=\begin{cases}-\frac{N\left({J}^{2}+{{\Gamma}}^{2}\right)}{4J}{\tau }_{f}\quad \hfill \\ -\frac{N{\Gamma}}{2}{\tau }_{f}.\quad \hfill \end{cases}\end{equation} \tag{ 39 }$$

The top solution in equation (39) corresponds to the two small-Γ SPs, while the bottom solution corresponds to the large-Γ SP. The double degeneracy of the former SP amounts to a factor of 2 multiplying one of the exponentials in (38): this does not contribute to the thermodynamic limit. Noticing that, for fixed J,

$$\begin{equation}\underset{N\to \infty }{\mathrm{lim}}\enspace \frac{1}{N}\enspace \mathrm{log}\sum _{\text{s.p.}}{\mathrm{e}}^{-N{\bar{S}}_{\mathrm{P}}\left({\Gamma}\right)}=-{\mathrm{min}}_{{\Gamma}}{\bar{S}}_{\mathrm{P}}\left({\Gamma}\right),\end{equation} \tag{ 40 }$$

where we defined the intensive quantities ${\bar{S}}_{\mathrm{P}}={S}_{\mathrm{P}}/N$, we obtain

$$\begin{equation}{{\epsilon}}_{\mathrm{G}}\left({\Gamma}\right)\approx {\mathrm{min}}_{{\Gamma}}\left({\bar{S}}_{\mathrm{P}}/{\tau }_{f}\right)=-{\mathrm{max}}_{{\Gamma}}\left(\frac{{J}^{2}+{{\Gamma}}^{2}}{4J},\frac{{\Gamma}}{2}\right),\end{equation} \tag{ 41 }$$

where the first solution is only valid for Γ < J, as discussed. Consistently with the findings of section 2.3, the LO result (41) is equal to the result of the MF approximation. The energy density of the true, entangled ground state is then retrieved upon including higher order terms in the expansion of (37).

Having obtained the LO MF approximation, we now wish to compute corrections beyond MF. In the presence of multiple SPs, it is customary to integrate Gaussian fluctuations about each SP and add up the relative contributions [66]. Here, we assume that the expansions about different SPs can be separately carried out and added up also for corrections beyond Gaussian. Let us discuss the conditions under which this procedure may be justified. Consider an integral whose integrand has several SPs. The expansion about each SP can be seen as a way of grouping contributions together: by expanding to higher and higher order, one progressively includes trajectories further and further from each SP. Adding up separate expansions around different SPs is then justified provided that there is no 'overlap': the trajectories included in one expansion are not significantly contributing to any of the others. A toy example showing this is provided in appendix C. In the present case, the requirement that there is no overlap is indeed satisfied: in the thermodynamic limit, equation (40) implies that only one expansion contributes for each value of Γ, and no double-counting can occur. Since the present discussion is not based on any model-specific assumption, the parameter Γ can here be understood as a general set of parameters specifying a given Hamiltonian. Additionally, in order to obtain finite results, each expansion should only be considered in the region of parameter space where it is convergent. This requirement can be physically understood as accounting for the breakdown of, for example, a large-coupling expansion in the small-coupling regime. With these caveats, let us carry out the full expansions as discussed; one has

$$\begin{equation}{{\epsilon}}_{\mathrm{G}}=\underset{{\tau }_{f}\to \infty }{\mathrm{lim}}\enspace \underset{N\to \infty }{\mathrm{lim}}\enspace \frac{1}{{\tau }_{f}N}\enspace \mathrm{log}\enspace \sum _{\text{s.p.}}{\mathrm{e}}^{-N{\bar{S}}_{\mathrm{P}}^{\prime }}\end{equation} \tag{ 42 }$$

where, for each SP, the quantity ${\bar{S}}_{\mathrm{P}}^{\prime }$ includes all contributions from higher order terms. Since the ground state energy density is finite and intensive, we expect ${\bar{S}}_{\mathrm{P}}^{\prime }$ to be independent of the system size. By equation (40), this means that

$$\begin{equation}{{\epsilon}}_{\mathrm{G}}={\mathrm{min}}_{{\Gamma}}\left({\bar{S}}_{\mathrm{P}}^{\prime }/{\tau }_{f}\right).\end{equation} \tag{ 43 }$$

Thus, the ground state energy is given by whichever of the series ${\bar{S}}_{\mathrm{P}}^{\prime }$, obtained by expanding around the different SPs, gives the lowest value of _G for a given value of the physical parameters. For example, for a quantum Ising model with fixed J, equation (43) implies that only one of the summands in (42) contributes for each choice of Γ. The above structure then suggests an interpretation of QPTs in terms of an abrupt change in which of the series ${\bar{S}}_{\mathrm{P}}^{\prime }$ yields the dominant contribution to the ground state energy. Let us denote by _±(Γ) the expansions corresponding to the two smallest ${\bar{S}}_{\mathrm{P}}^{\prime }$ in a given range of the parameter Γ. In the thermodynamic limit, due to the minimum function in (43), the functional form of _G changes abruptly at the value Γ_c such that ₋ < ₊ if Γ < Γ_c and ₋ > ₊ if Γ > Γ_c, so that the value Γ_c can be identified as the quantum critical point. We refer to this condition as crossing of the two series, although both may diverge at the critical point Γ_c itself; see appendix C. Equation (42) thus implies that the GS energy density can only be non-analytic in the thermodynamic limit, and only at the point where two different series cross in the above-defined sense. This expected result is here retrieved purely on the basis of the disentanglement formalism, where it emerges naturally as a consequence of equations (36), (40) and (42). In the next section, we will show that for the Ising model the expansions around the small-Γ and large-Γ SPs give rise to series in Γ/J and J/Γ respectively. These can be identified as the small-Γ and large-Γ perturbative expansions of the ground state energy. To benchmark our analytical approach, we shall exploit the exact solvability of the Ising chain, whose ground state energy density is given by [65]

$$\begin{equation}{{\epsilon}}_{\mathrm{G}}=-\frac{2{\Gamma}+J}{2\pi }E\left(\frac{8J{\Gamma}}{{\left(J+2{\Gamma}\right)}^{2}}\right),\end{equation} \tag{ 44 }$$

in terms of the complete elliptic integral of the second kind E. Expanding equation (44) for small or large Γ, one finds

$$\begin{equation}{{\epsilon}}_{\mathrm{G}}=\begin{cases}-\frac{J}{4}-\frac{{{\Gamma}}^{2}}{4J}-\frac{{{\Gamma}}^{4}}{16{J}^{3}}-\frac{{{\Gamma}}^{6}}{16{J}^{5}}-\dots ,\quad \hfill \\ -\frac{{\Gamma}}{2}-\frac{{J}^{2}}{32{\Gamma}}-\frac{{J}^{4}}{2048{{\Gamma}}^{3}}-\frac{{J}^{6}}{32\enspace 768{{\Gamma}}^{5}}-\dots ,\quad \hfill \end{cases}\end{equation} \tag{ 45 }$$

respectively. When all terms are resummed, the above perturbative series do indeed cross only at the critical point Γ_c = J/2; see appendix C. The small-Γ expansion is seen to be divergent for Γ > Γ_c; therefore, the relative terms will not contribute in this regime. Similarly, the large-Γ series does not contribute when Γ < Γ_c. Below, we will show how small- and large-Γ expansions for _G are obtained from the disentanglement approach.

The evaluation of corrections to the GS energy beyond LO according to equation (42) is carried out by expanding the functional integral representation of (37) around the SPs and applying Wick's theorem. This gives rise to a set of Feynman rules, which we illustrate in this section by considering the quantum Ising chain. Since the expansions around the different SPs take the same functional form, here we work in full generality, specializing only the final results to each SP. In order to reveal the dependence of each term on the physical parameters Γ and J, it is convenient to introduce dimensionless times $\bar{\tau }={\gamma }_{\mathrm{P}}\tau$, where γ_P will eventually be set to the appropriate plateau value (33) for each SP. The noise action becomes

$$\begin{equation}{S}_{0}\equiv \frac{J}{4{\gamma }_{\mathrm{P}}}{\int }_{0}^{{\bar{\tau }}_{f}}\mathrm{d}\tau \sum _{ij}{\left({\mathcal{J}}^{-1}\right)}_{ij}{\varphi }_{i}\left(\tau \right){\varphi }_{j}\left(\tau \right)\end{equation} \tag{ 46 }$$

and the effective Loschmidt action is given by

$$\begin{equation}S={S}_{0}-{\int }_{0}^{{\bar{\tau }}_{f}}\mathrm{d}\bar{\tau }\left[\frac{{\Gamma}}{2{\gamma }_{\mathrm{P}}}\sum _{i}{\xi }_{i}^{+}\left(\bar{\tau }\right)+\frac{J}{{\gamma }_{\mathrm{P}}}\sum _{i}{\varphi }_{i}\left(\bar{\tau }\right)\right].\end{equation} \tag{ 47 }$$

It is convenient to separately consider the variations of (47) which involve functional derivatives of ${\xi }_{i}^{+}$ and the term originating from the noise action S₀: as we shall see, the latter provides a simple and physically appealing propagator ² . Thus, for each SP we expand the action (47) as

$$\begin{equation}S={S}_{\mathrm{P}}+{S}_{0}-\sum _{n=2}^{\infty }{T}^{\left(n\right)}\end{equation} \tag{ 48 }$$

where we defined

$$\begin{equation}{T}^{\left(n\right)}\equiv \frac{1}{n!}\sum _{i}{\int }_{0}^{{\bar{\tau }}_{f}}\dots {\int }_{0}^{{\bar{\tau }}_{f}}{S}^{\left(n\right)}\left({\bar{\tau }}_{1},\dots ,{\bar{\tau }}_{n}\right){\varphi }_{i}\left({\bar{\tau }}_{1}\right)\dots {\varphi }_{i}\left({\bar{\tau }}_{n}\right)\mathrm{d}{\bar{\tau }}_{1}\dots \mathrm{d}{\bar{\tau }}_{n},\end{equation} \tag{ 49 }$$

$$\begin{equation}{S}^{\left(n\right)}\left({\bar{\tau }}_{1},\dots ,{\bar{\tau }}_{n}\right)\equiv -\frac{{\Gamma}}{2{\gamma }_{\mathrm{P}}}{\int }_{0}^{{\bar{\tau }}_{f}}\mathrm{d}\bar{\tau }{\left.\frac{{\delta }^{n}{\xi }_{i}^{+}\left(\bar{\tau }\right)}{\delta {\varphi }_{i}\left({\bar{\tau }}_{1}\right)\dots \delta {\varphi }_{i}\left({\bar{\tau }}_{n}\right)}\right\vert }_{\mathrm{P}},\end{equation} \tag{ 50 }$$

and we exploited translational invariance and $\delta {\xi }_{i}^{+}/\delta {\varphi }_{j}\propto {\delta }_{ij}$. The functional integral (24) can be expanded as

$$\begin{equation}A\left({\tau }_{f}\right)=\sum _{\text{s.p.}}{\mathrm{e}}^{-{S}_{\mathrm{P}}}{\left\langle \sum _{m=0}^{\infty }\frac{{\left({\sum }_{n=2}^{\infty }\;{T}^{\left(n\right)}\right)}^{m}}{m!}\right\rangle }_{0},\end{equation} \tag{ 51 }$$

where the notation ⟨...⟩₀ denotes averaging with respect to the noise action (46). Each term in equation (51) can be evaluated using Wick's theorem. The propagator Δ, which accounts for interactions in the system, can be read off from the quadratic action (46) and is found to be proportional to the interaction matrix:

$$\begin{equation}{{\Delta}}_{ij}\left(\bar{\tau },{\bar{\tau }}^{\prime }\right)=2\frac{{\gamma }_{\mathrm{P}}}{J}{\mathcal{J}}_{ij}\delta \left(\bar{\tau }-{\bar{\tau }}^{\prime }\right).\end{equation} \tag{ 52 }$$

The series obtained from (51) can then be formally re-exponentiated, giving equation (42). Consider the averages ⟨...⟩₀ in equation (51). We define the order of a term ${\langle {T}^{\left({n}_{1}\right)}\dots {T}^{\left({n}_{m}\right)}\rangle }_{0}$ to be $l={\sum }_{j=1}^{m}{n}_{j}$. Wick's theorem implies that terms of odd order vanish identically, while terms of even order are obtained by summing over all the possible replacements of pairs of fields ${\phi }_{i}\left({\bar{\tau }}_{i}\right)$, ${\phi }_{j}\left({\bar{\tau }}_{j}\right)$ by propagators ${{\Delta}}_{ij}\left({\bar{\tau }}_{i},{\bar{\tau }}_{j}\right)$. The evaluation of a given term in (51) is simplified by means of a diagrammatic representation and the associated Feynman rules:

• Each T⁽ⁿ⁾ provides a vertex and contributes a factor S⁽ⁿ⁾/n!. Diagrammatically, a vertex is represented as n points arranged inside a box. Each vertex is labeled by a site index, e.g. j. Individual points belonging to a given vertex are additionally distinguished by a unique time label, ${\bar{\tau }}_{{j}_{1}}\dots {\bar{\tau }}_{{j}_{n}}$.
• One then sums over all possible ways of joining pairs of points by lines; a line joining the points labeled by $\left(j,{\bar{\tau }}_{j}\right)$, and $\left(k,{\bar{\tau }}_{k}\right)$ gives a propagator ${{\Delta}}_{jk}\left({\bar{\tau }}_{j},{\bar{\tau }}_{k}\right)$.
• The resulting quantity is then integrated over all times ${\bar{\tau }}_{i}$; the integrals run between 0 and ${\bar{\tau }}_{f}$.
• Finally, all site indices are summed over.

Examples of this diagrammatic representation are given in figures 1–3, discussed below. In more usual field theories, such as ϕ⁴, vertices are typically represented by single points [67]; in the above rules, this would correspond to setting all ${\bar{\tau }}_{{j}_{i}}$ to the same value. The fact that here vertices consist of separate points is due to the non-locality in time of the action (48). To simplify the evaluation of higher order terms, we identify two classes of diagrams which do not contribute. The first class includes diagrams where lines join points within the same vertex: these are self-interaction diagrams. An example, originating from ${\langle {T}^{\left(3\right)}{T}^{\left(3\right)}{T}^{\left(2\right)}\rangle }_{0}$, is shown in figure 1(a). Such diagrams feature at least one term of the form ${{\Delta}}_{ii}\propto {\mathcal{J}}_{ii}$, which is identically zero for the Hamiltonian (20) at hand. The second class of non-contributing terms includes disconnected diagrams, in which the vertices joined by internal lines form disjointed clusters. This class includes the diagram in figure 1(b), which is produced by ${\langle {T}^{\left(4\right)}{T}^{\left(4\right)}\rangle }_{0}$. It is easy to see that, due to the Feynman rules and the form of the propagator (52), a connected cluster of vertices gives a contribution proportional to the system size N. A term with m disconnected clusters of vertices is then proportional to N^m . The origin of these terms can be understood by considering the expansion

$$\begin{equation}{\mathrm{e}}^{-N{\bar{S}}^{\prime }}=1-N{\bar{S}}^{\prime }+\frac{1}{2!}{N}^{2}{\left({\bar{S}}^{\prime }\right)}^{2}+\dots .\end{equation} \tag{ 53 }$$

The terms ⟨...⟩₀ in (51), obtained from Wick's theorem, correspond to the right-hand side of (53). The ground state energy must be intensive and finite in the thermodynamic limit; this implies that terms proportional to higher powers of N must cancel out when exponentiating the series in (51) to obtain equation (42). This is precisely what happens in equation (53). Since the desideratum here is ${\bar{S}}^{\prime }$, we only need to consider terms proportional to N: these are given by connected diagrams. Finally, we note that for the quantum Ising model all diagrams with an odd number m of vertices T⁽²⁾ vanish: any such diagram is either self-interacting (figure 1(c)), or it gives rise to a term $\propto \mathrm{T}\mathrm{r}\enspace {\mathcal{J}}^{m}$, which vanishes for ${\mathcal{J}}_{ij}=\left({\delta }_{i,j+1}+{\delta }_{i+1,j}\right)/2$ (figure 1(d)).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In section 3.5, we apply the Feynman rules derived in this section to compute higher order corrections to the ground state energy of the quantum Ising chain. Before turning to this explicit example, in the next section we complete our theoretical overview by considering how the terms in equation (51) depend on the physical parameters of the model, elucidating the relation between the expansion about the SPs and perturbation theory.

In order to understand the nature of the terms produced by the expansion (51) we need to consider the higher variations of the action, S⁽ⁿ⁾ with n ⩾ 2. It is convenient to compute these variations by initially imposing a time ordering ${\bar{\tau }}_{n}{ >}\cdots { >}{\bar{\tau }}_{1}$, and then symmetrizing the result with respect to the times ${\bar{\tau }}_{i}$. With said ordering, one obtains for the second variation

$$\begin{align}\hfill {S}^{\left(2\right)}\left({\bar{\tau }}_{1},{\bar{\tau }}_{2}\right)& =-{\left(\frac{J}{{\gamma }_{\mathrm{P}}}\right)}^{2}{\int }_{0}^{{\bar{\tau }}_{f}}\frac{{\Gamma}}{{\gamma }_{\mathrm{P}}}{\bar{{\Xi}}}_{i}\left({\bar{s}}_{1},{\bar{\tau }}_{1}\right)\hfill \\ \hfill & \quad {\times}{\left.\left[\theta \left({\bar{s}}_{1}-{\bar{\tau }}_{2}\right)-\frac{{\Gamma}}{{\gamma }_{\mathrm{P}}}{\int }_{{\bar{\tau }}_{2}}^{{\bar{s}}_{1}}{\bar{{\Xi}}}_{i}\left({\bar{s}}_{2},{\bar{\tau }}_{2}\right)\mathrm{d}{\bar{s}}_{2}\right]\mathrm{d}{\bar{s}}_{1}\right\vert }_{\mathrm{P}},\hfill \end{align} \tag{ 54 }$$

where we defined ${\Xi}\left(\bar{s},{\bar{\tau }}_{1}\right)\equiv \frac{J}{{\gamma }_{\mathrm{P}}}\bar{{\Xi}}\left(\bar{s},{\bar{\tau }}_{1}\right)$ to make the dependence on physical parameters manifest. Equation (54) shows that all variations S⁽ⁿ⁾ with n ⩾ 2 can be expressed in terms of integrals of the first variation Ξ_i . Schematically, one obtains S⁽ⁿ⁺¹⁾ from S⁽ⁿ⁾ by summing over all possible ways of replacing

$$\begin{equation*}\bar{{\Xi}}\to \theta -\frac{{\Gamma}}{{\gamma }_{\mathrm{P}}}\int \bar{{\Xi}}\end{equation*}$$

and multiplying by J/γ_P. When evaluated at the SP, each $\bar{{\Xi}}$ gives a factor of ξ_P and an exponential depending on the dimensionless times ${\bar{\tau }}_{i}$ only. Therefore, the nth variation (with n ⩾ 2) evaluated at the plateau must be of the form

$$\begin{equation}{S}^{\left(n\right)}\left({\bar{\tau }}_{1},\dots ,{\bar{\tau }}_{n}\right)={\left(\frac{J}{{\gamma }_{\mathrm{P}}}\right)}^{n}\sum _{m=1}^{n}{C}_{n,m}\left(\bar{\tau }\right){\left(\frac{{\Gamma}{\xi }_{\mathrm{P}}^{+}}{{\gamma }_{\mathrm{P}}}\right)}^{m},\end{equation} \tag{ 55 }$$

where ${C}_{n,m}\left(\bar{\tau }\right)$ are dimensionless functions depending only on the ${\bar{\tau }}_{i}$, which do not involve any factor of Γ, J or γ_P. Using the form (55) of higher variations, it is possible to determine the structure of the terms in equation (51). Terms of odd order l = 2m + 1 do not contribute due to Wick's theorem; see the discussion in section 3.3. Any given term of even order l = 2m features a product of variations, whose orders add up to 2m, and m propagators, each of which carries a factor γ_P/J. Bringing everything together and substituting the SP values (33), a general 2mth order term ${\mathcal{T}}^{2m}$ can be written as

$$\begin{equation}{\mathcal{T}}^{2m}=\begin{cases}\sum _{n=1}^{\infty }{\bar{C}}_{n,2m}{\left(\frac{{\Gamma}}{J}\right)}^{2n}\quad \hfill \\ {\bar{D}}_{2m}{\left(\frac{J}{{\Gamma}}\right)}^{m}\quad \hfill \end{cases}\end{equation} \tag{ 56 }$$

where the top and bottom solutions refer to the small-Γ and large-Γ expansion respectively, and ${\bar{C}}_{n,2m}$ and ${\bar{D}}_{2m}$ are dimensionless constants which do not depend on J or Γ. The series of even powers of (Γ/J) in the former case arises from Taylor expanding ${\Gamma}{\xi }_{\mathrm{P}}^{+}/{\gamma }_{\mathrm{P}}$. Equation (56) thus shows that the expansions around the SPs give rise to series in (Γ/J)² and J/Γ. These can be identified with the perturbative series as follows. Equation (42) is valid for any value of Γ and, due to the thermodynamic relation (43), only one expansion at a time contributes. Consider the Γ → 0 limit; in this case, the LO term (39) of the small-Γ expansion gives the exact value of the ground state energy. For finite but sufficiently small Γ/J ≪ 1, the small-Γ expansion will still be the dominant one and give the ground state energy _G. The small-Γ expansion must therefore be equal to the Γ/J perturbative series for _G, as they are both series in Γ/J and they both add up to _G. A symmetric argument holds for the large-Γ series in the corresponding limit. Although we illustrated this result for the Ising model, the correspondence is expected be more generally valid: both expansions yield the exact free energy (Γ) in the strong-coupling (Γ = 0) or strong-field (Γ → ∞) limits, and by smoothness of (Γ) it can be expected that these expansions will continue to give the same result provided no phase transition is crossed. In the case of the Ising model, the argument applies either side of the transition, as discussed above.

Equation (56) further shows that the large-Γ expansion (bottom case) is in order-by-order correspondence to the perturbative series: terms of order 2m are proportional to (Γ/J)^m . On the other hand, the small-Γ expansion (top case of equation (56)) is not in one-to-one correspondence to perturbation theory: one needs to sum over m in order to retrieve the perturbative series in Γ/J, since each of the terms in equation (56) may in principle contain various powers of Γ/J. Such different behavior of the two expansions is due to the nature of the plateau configuration or, equivalently, the MF ground state. For large Γ, this is just the Γ = ∞ ground state |⇒⟩; the large-Γ expansion is thus equivalent order-by-order to the perturbative series around Γ = ∞. On the other hand, for small Γ the MF ground state is not simply given by the Γ = 0 ground state |⇓⟩. Consider for instance the MF magnetization, given in appendix D; this can be expanded as a Taylor series featuring all even powers of (Γ/J). An expansion around the MF ground state is therefore not expected to be in order-to-order correspondence with a perturbative expansion around Γ = 0.

One more comment is due concerning even and odd powers in the two expansions. From expanding the exact ground state energy of the quantum Ising chain as in (45), we see that the perturbative expression for _G/J around Γ = 0 features only even powers of Γ/J and, similarly, the perturbative expansion of _G/Γ around Γ = ∞ contains only even powers of J/Γ. This result is immediately retrieved from equation (56) for the small-Γ expansion, and is due to spontaneous symmetry breaking at the MF level. However, odd powers of J/Γ are not excluded a priori in the large-Γ expansion. The necessary cancellation must therefore originate from the vanishing of the ${\bar{D}}_{2m}$ coefficient in (56) when m is odd. We explicitly show an example of such cancellation when computing higher-order corrections to _G in section 3.5.

We have thus determined the structure of the terms produced by expanding about the SPs, and clarified the relation of such expansions to perturbation theory. In summary, the small-Γ and large-Γ expansions, taken as a whole, are respectively equal to the full perturbative series around Γ = 0 and Γ = ∞. This correspondence is satisfied order-by-order for the large-Γ expansion, but only when resumming the whole series for the small-Γ expansion. This completes an overview of the field theoretical approach; in the next section, we apply the concepts discussed so far to compute corrections to the ground state energy of the quantum Ising chain.

In order to illustrate the machinery introduced in the previous sections, we compute the next-to-leading order (NLO) and next-to-next-to-leading order (NNLO) corrections to the ground state energy density of the quantum Ising chain. The lowest order correction is naively given by ${\langle {T}^{\left(2\right)}\rangle }_{0}$; this term however vanishes, since the corresponding diagram is self-interacting: see the discussion in section 3.3. The NLO correction is then given by the next-higher term, which is of order four:

$$\begin{equation}{\mathcal{T}}^{\left(4\right)}=\frac{1}{2}{\langle {T}^{\left(2\right)}{T}^{\left(2\right)}\rangle }_{0}+{\langle {T}^{\left(4\right)}\rangle }_{0}.\end{equation} \tag{ 57 }$$

The second term on the right-hand side of equation (57) vanishes similarly to ${\langle {T}^{\left(2\right)}\rangle }_{0}$; the remaining term corresponds to the diagrams in figure 2, and gives

$$\begin{equation}\frac{1}{2}{\langle {T}^{\left(2\right)}{T}^{\left(2\right)}\rangle }_{0}={\bar{\tau }}_{f}N\begin{cases}\frac{{{\Gamma}}^{4}}{32{J}^{4}}\quad \hfill \\ \frac{{J}^{2}}{32{{\Gamma}}^{2}}\quad \hfill \end{cases}\end{equation} \tag{ 58 }$$

where the top and bottom results are obtained from the small-Γ and large-Γ SPs respectively. Thus, the NLO approximation to the ground state energy to is given by

$$\begin{equation}{{\epsilon}}_{\mathrm{G}}\approx -{\mathrm{max}}_{{\Gamma}}\left(\frac{J}{4}+\frac{{{\Gamma}}^{2}}{4J}+\frac{{{\Gamma}}^{4}}{32{J}^{3}},\frac{{\Gamma}}{2}+\frac{{J}^{2}}{32{\Gamma}}\right).\end{equation} \tag{ 59 }$$

As discussed in section 3.2, each of the two series in equation (59) can only be considered within its radius of convergence. In practice, one typically does not have access to the full series; the radius of convergence can then be estimated by imposing that each term be smaller than the lower-order one. In the present case, this criterion indicates that the small-Γ series is valid for Γ < 1, while the large-Γ series is valid for Γ > 1/4. Including the NLO correction as discussed provides an improvement over the LO approximation for all values of Γ; see figure 4. Consistently with the discussion in section 3.4, equation (59) matches the result of second-order perturbation theory about Γ = ∞, which is in one-to-one correspondence with the expansion around the large-Γ SP. On the other hand, in order to match the term of order Γ⁴/J³ of the perturbative series in (45), one needs to include higher order contributions from the small-Γ SP expansion: this is again consistent with the earlier discussion. The next-higher correction to the ground state energy, NNLO, is of order 6, and is given by

$$\begin{equation}{\mathcal{T}}^{\left(6\right)}=\frac{1}{3!}{\langle {T}^{\left(2\right)}{T}^{\left(2\right)}{T}^{\left(2\right)}\rangle }_{0}+\frac{1}{2!}{\langle {T}^{\left(3\right)}{T}^{\left(3\right)}\rangle }_{0}.\end{equation} \tag{ 60 }$$

The first term in equation (60) vanishes because it features an odd number of T⁽²⁾ vertices; see the discussion in section 3.3 and in particular figures 1(c) and (d). The non-vanishing diagrams are shown in figure 3. They can be evaluated to give

$$\begin{equation}{\mathcal{T}}^{\left(6\right)}={\bar{\tau }}_{f}N\begin{cases}\frac{{{\Gamma}}^{4}}{64{J}^{4}}-\frac{{{\Gamma}}^{6}}{64{J}^{6}}\quad \hfill \\ 0\quad \hfill \end{cases}\end{equation} \tag{ 61 }$$

where again the top result corresponds to the small-Γ expansion and the bottom result to the large-Γ expansion. Equation (61) shows that the NNLO correction from the large-Γ expansion vanishes. This was anticipated in section 3.4, and is due to the fact that the SP expansions and the perturbative series must coincide order-by-order; by equation (56), the large-Γ NNLO correction would be proportional to J³, but no such term appears in the perturbative series (45): the coefficient multiplying J³ must therefore vanish. Including the NNLO term (61) leads to a further improvement in the approximation to the GS energy, as shown in the inset of figure 4.

Zoom In Zoom Out Reset image size

We have thus illustrated how the disentanglement formalism provides a method to analytically approximate ground state expectation values. The LO SP result corresponds to MF, while analytical corrections beyond MF can be systematically obtained by means of Feynman diagrams, which we explicitly showed. For simplicity, in this section we focused on the ground state energy, but an analogous procedure can be carried out for other observables by expanding the appropriate effective action ${S}_{\mathcal{O}}$, defined as in equation (18). The machinery introduced above is qualitatively different from other existing diagrammatic techniques. Diagrammatic expansions for spin systems conventionally make use of a formal mapping to an auxiliary fermionic system with imaginary chemical potential [68]; see e.g. [69–71] for recent applications. In contrast, the present field theory is directly formulated in spin language. Furthermore, the present method does not yield a perturbative expansion about a classical or non-interacting limit, respectively given by Γ = 0 and J = 0 for the Ising model, but an expansion about the MF solution. These are in general different, e.g. the MF magnetization for the D-dimensional Ising model is given by $\sqrt{{D}^{2}{J}^{2}-{{\Gamma}}^{2}}/2$ for Γ < JD, whose small −Γ expansion features all even powers of Γ; this cannot be matched by any order of perturbation theory.

Beside the analytical field theoretical formalism outlined in this section, the disentanglement method can alternatively be used as a numerical tool; we discuss this approach in the following section.

Knowledge of the SP trajectory can be used to implement an importance sampling numerical algorithm. For this application, it is convenient to work with the diagonal form (11) of the noise action, involving the fields ${\phi }_{i}^{a}$. The corresponding SP values can be readily determined using ${\phi }_{i}^{a}{\vert }_{\text{SP}}={\sum }_{bj}\;{O}_{ij}^{ab}{\varphi }_{j}^{b}{\vert }_{\text{SP}}$. The key step of the proposed approach consists in using ${\phi }_{\text{SP}}\equiv \left\{{\phi }_{i}^{a}{\vert }_{\text{SP}}\right\}$ to perform the change of variables

$$\begin{equation}{\phi }_{i}^{a}={\left({\phi }_{\text{SP}}\right)}_{i}^{a}+{\phi }_{i}^{\prime a}\end{equation} \tag{ 62 }$$

in the functional integral for a given observable:

$$\begin{equation}\langle \hat{\mathcal{O}}\rangle ={\mathrm{e}}^{-{S}_{0}\left[{\phi }_{\text{SP}}\right]}\int \mathcal{D}{\phi }^{\prime }\enspace {\mathrm{e}}^{-{S}_{0}\left[{\phi }^{\prime }\right]}\enspace {\mathrm{e}}^{-\int \mathrm{d}\tau {\phi }^{\prime }\left(\tau \right)\cdot {\phi }_{\text{SP}}\left(\tau \right)}{f}_{\mathcal{O}}\left[{\phi }_{\text{SP}}+{\phi }^{\prime }\right],\end{equation} \tag{ 63 }$$

where ${\phi }^{\prime }\equiv \left\{{\phi }_{i}^{\prime a}\right\}$ and ${\phi }_{\text{SP}}\cdot {\phi }^{\prime }\equiv {\sum }_{ia}{\left({\phi }_{\text{SP}}\right)}_{i}^{a}{\phi }_{i}^{a}$. Due to the Gaussianity of the noise action S₀, equation (63) can be evaluated numerically in the spirit of the stochastic approach of [58, 59]; this amounts to averaging a biased function over realizations of Gaussian-distributed stochastic processes ϕ':

$$\begin{equation}\langle \hat{\mathcal{O}}\rangle ={\mathrm{e}}^{-{S}_{0}\left[{\phi }_{\text{SP}}\right]}{\langle {\mathrm{e}}^{-\int \mathrm{d}\tau {\phi }^{\prime }\left(\tau \right)\cdot {\phi }_{\text{SP}}\left(\tau \right)}{f}_{\mathcal{O}}\left[{\phi }_{\text{SP}}+{\phi }^{\prime }\right]\rangle }_{{\phi }^{\prime }}.\end{equation} \tag{ 64 }$$

By construction, the stochastic processes ϕ' featured in (64) are fluctuations about the dominant SP trajectories. In contrast, when equation (18) is sampled directly according to (11), trajectories close to ϕ(τ) = 0 are sampled preferentially, even though they may give a small contribution to the integral. Compared to more usual path integral approaches, here we do not truncate to a given order in the fluctuations, so that equation (63) does not constitute a semiclassical approximation. Instead, a change of variables is used to bias the sampling toward important trajectories; the exactness of equation (18) is thus fully preserved in equation (63). This change of variables can be seen as a particular measure (or Girsanov) transformation [57, 72]; in the context of stochastic processes, this constitutes the continuum version of importance sampling.

To illustrate our method, we apply the measure transformation approach to the D-dimensional quantum Ising model (20) for D ∈ {1, 2, 3} by numerically computing different observables from stochastic simulations. In our numerical simulations and in the remainder of this section we set J = 1. In the stochastic approach, ground state expectation values are computed according to equation (4). By appropriately choosing ${\hat{U}}_{0}$ in equation (3), any initial state |ψ₀⟩ can be considered. Following the discussion of section 2.3, it is convenient to choose the initial state to be the mean-field ground state for the desired value of Γ, |MF⟩. As anticipated, this is equivalent to initializing the system at the plateau SP configuration, ${\xi }_{i}^{+}\left(0\right)={\xi }_{\mathrm{P}}^{+}$, ${\xi }_{i}^{z}\left(0\right)=\mathrm{log}\left(1+\vert {\xi }_{\mathrm{P}}^{+}{\vert }^{2}\right)$. For observables computed from (4), the plateau values are fixed points of the SP equation, such that one may perform the change of variables (62) with ϕ_SP(τ) = ϕ_P; see appendix C. The subsequent imaginary time evolution then projects the wavefunction from the MF to the true ground state. The SDEs implementing the imaginary time evolution are solved using the Euler scheme [57]. The statistical uncertainty on each quantity obtained as a classical average is computed by partitioning the data set into n_B batches of independent simulations. We use n_B = 100 unless otherwise stated. The data within each batch are averaged, and the standard deviation σ of each quantity over the n_B batches is then computed. The uncertainty on the mean is estimated as the standard error $\sigma /\sqrt{{n}_{\mathrm{B}}}$. Since observables are obtained from the ratio (4), we apply the appropriate uncertainty propagation formula

$$\begin{equation}\frac{{\sigma }^{2}\left(X/Y\right)}{{\left(X/Y\right)}^{2}}=\frac{{\sigma }^{2}\left(X\right)}{{X}^{2}}+\frac{{\sigma }^{2}\left(Y\right)}{{Y}^{2}}-2\frac{\mathrm{c}\mathrm{o}\mathrm{v}\left(X,Y\right)}{XY},\end{equation} \tag{ 65 }$$

where cov(X, Y) is the covariance of X, Y. It is also convenient to exploit the real-valuedness of the numerator and denominator of (4) to consider only the real parts of quantities obtained from averaging. We estimate the error on a given observable $\mathcal{O}$ to be $\sigma \left(\mathcal{O}\right)$ computed as above.

We begin by considering the D = 1 case, corresponding to the quantum Ising chain. We consider the imaginary time evolution of the ground state longitudinal magnetization ${\mathcal{M}}_{z}$, transverse magnetization ${\mathcal{M}}_{x}\equiv {\sum }_{i=1}^{N}\langle {\hat{S}}_{i}^{x}\rangle /N$ and nearest-neighbor longitudinal correlations ${C}_{zz}\equiv {\sum }_{\langle ij\rangle }\langle {\hat{S}}_{i}^{z}{\hat{S}}_{j}^{z}\rangle /N$. As per our general discussion, these quantities are given by equation (13), where the numerator includes the appropriate stochastic function for each observable. The stochastic functions are given by equation (16) for ${\mathcal{M}}_{z}$ and by

$$\begin{equation}{F}_{{\mathcal{M}}_{x}}=\frac{{F}_{\mathbb{1}}}{N}\sum _{i}\frac{{\xi }_{\mathrm{f},i}^{+}\left(\tau \right)+{\xi }_{\mathrm{b},i}^{+{\ast}}\left(\tau \right)}{1+{\xi }_{\mathrm{f},i}^{+}\left(\tau \right){\xi }_{\mathrm{b},i}^{+{\ast}}\left(\tau \right)},\end{equation} \tag{ 66 }$$

$$\begin{equation}{F}_{{C}_{zz}}=\frac{{F}_{\mathbb{1}}}{N}\sum _{\langle ij\rangle }\left(\frac{1-{\xi }_{\mathrm{f},i}^{+}\left(\tau \right){\xi }_{\mathrm{b},i}^{+{\ast}}\left(\tau \right)}{1+{\xi }_{\mathrm{f},i}^{+}\left(\tau \right){\xi }_{\mathrm{b},i}^{+{\ast}}\left(\tau \right)}\right)\left(\frac{1-{\xi }_{\mathrm{f},j}^{+}\left(\tau \right){\xi }_{\mathrm{b},j}^{+{\ast}}\left(\tau \right)}{1+{\xi }_{\mathrm{f},j}^{+}\left(\tau \right){\xi }_{\mathrm{b},j}^{+{\ast}}\left(\tau \right)}\right),\end{equation} \tag{ 67 }$$

for ${\mathcal{M}}_{x}$ and C_zz respectively [59]. The ground state energy density is obtained from equations (66) and (67) as

$$\begin{equation}{{\epsilon}}_{\mathrm{G}}=-{\Gamma}{\mathcal{M}}_{x}-J{C}_{zz}.\end{equation} \tag{ 68 }$$

In figure 5 we compare our numerical results for a system of size N = 101 to imaginary time evolution performed directly in the thermodynamic limit using iTEBD [48]. We find excellent agreement across the imaginary time range we consider. The error on the energy estimate (τ_f ) we obtain for τ_f = 4 is approximately 1.5 × 10⁻⁵ relative to the exact ground state result. To show the improvement of importance sampling according to equation (63) over direct sampling using the naive measure (11), in figure 6 we compare the performance of the two approaches. We fix the physical parameters to N = 15, Γ = 0.4 and compute the Euclidean time evolution of the energy using the same time step and number of simulations; the results obtained from direct and importance sampling are shown in panels (a) and (b) respectively. It is clear that the importance sampling algorithm produces far better results for the same computational cost. This is further discussed in section 4.3, where we study the behavior of fluctuations. The importance sampling algorithm can be equally applied to higher dimensional systems, for which analytical solutions are typically not available. The relevant stochastic formulae take a similar form to the one dimensional case [58]. For instance, the stochastic function for the normalization and the longitudinal magnetization are given by equations (15) and (16) respectively, where the sums and products are performed over all lattice sites. Similarly, the measure transformation is carried out in complete analogy to the one-dimensional case; see appendix C for further details. In figure 7, we consider the 2D quantum Ising model, comparing the results obtained from importance sampling and from exact diagonalization (ED) performed with the QuSpin package [73]. We consider a 5 × 5 system for Γ = 1, performing Euclidean time evolution from the MF ground state. Again, we find excellent agreement between our result and ED; the relative error on our estimate for the ground state energy is within 3 × 10⁻⁵. Finally, in figure 8 we consider the quantum Ising model in three spatial dimensions. Our results are again in good agreement with ED for a system of size 3 × 3 × 3; for the chosen stopping time τ_f , we obtain a relative error on the GS energy within 10⁻³. We observe that the stopping time τ_f that can be accessed for a given number of simulations decreases with the dimensionality of the system, due to the faster growth of fluctuations; this is due to the greater connectivity in higher dimensions, and is further investigated in the next section.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Having demonstrated the applicability of the method to higher dimensional systems, we now turn to investigating its numerical performance, quantitatively comparing the direct and importance sampling schemes. Fluctuations in the stochastic quantities play a central role: for a given number of simulations, the growth of fluctuations ultimately determines the time scale beyond which physical results are not correctly reproduced. Therefore, an increasing number of simulations is needed as the stopping time is increased. The central limit theorem implies that the fluctuations in the observable $\mathcal{O}$ computed from the stochastic approach are determined by the variance σ² of the corresponding stochastic quantity ${f}_{\mathcal{O}}$ [59]:

$$\begin{equation}{\sigma }^{2}\left({f}_{\mathcal{O}}\right)\equiv {\langle \vert {f}_{\mathcal{O}}{\vert }^{2}\rangle }_{\phi }-\vert {\langle {f}_{\mathcal{O}}\rangle }_{\phi }{\vert }^{2}.\end{equation} \tag{ 69 }$$

The variance σ² is therefore directly related to the number of simulations required to obtain a given accuracy. We illustrate this by considering the normalization function (15): the behavior of this quantity is found to be representative of other observables, due to the similar functional form of the corresponding stochastic functions; see for example equation (16). In the classical case of the D-dimensional Ising model with Γ = 0, the SDEs (21) are exactly solvable. For direct sampling, one obtains

$$\begin{equation}{\sigma }^{2}\left(\tau \right)={\mathrm{e}}^{2NDJ\tau }-{\mathrm{e}}^{NDJ\tau }\end{equation} \tag{ 70 }$$

where ND is the total number of interactions in the system. In contrast, the variance σ² for the importance sampling scheme vanishes identically and a single trajectory is sufficient to give the exact result. For finite Γ, the behavior of fluctuations in the two approaches can be investigated numerically. As shown in figure 9, we find that this is captured by the functional form

$$\begin{equation}{\sigma }^{2}=\alpha {\text{e}}^{\beta \tau },\end{equation} \tag{ 71 }$$

with β ≈ 2DN. Thus, the exponential growth of fluctuations with N, D and τ, which we found for direct sampling in the classical case Γ = 0, survives for finite Γ, and also applies to importance sampling. This is consistent with the numerical analysis carried out in [59] for real time evolution, and with the argument of [24, 64] suggesting that a large deviation principle (LDP) may be at play with respect to the system size N. However, direct and importance sampling differ substantially in the prefactor α multiplying the exponential. We find that α depends heavily on Γ, as shown in the tables of figure 9. For direct sampling, one has α = O(1) for all Γ. On the other hand, for importance sampling, α gradually increases from zero as Γ is increased. For small to intermediate field strengths Γ ≈ O(1), α can be orders of magnitude smaller for importance sampling than for direct sampling. Thus, the importance sampling scheme can strongly suppress the growth of fluctuations with time and the system size. This significantly extends the regime of applicability of the stochastic method before fluctuations become sizable, although it does not eliminate their exponential growth. It would be interesting to clarify the relation between stochastic fluctuations and entanglement. The importance sampling method completely eliminates fluctuations when the true ground state is a product state and does not have entanglement, as in the classical limit. Thus, both the presence of residual fluctuations and the growth of entanglement signal the departure from a product state; whether a direct connection between these exists will be investigated in future work.

Zoom In Zoom Out Reset image size

The computational cost of the numerical stochastic approach is mainly determined by the growth of fluctuations. The runtime of a given stochastic simulation scales linearly with N, τ_f and inversely with the time step Δτ_f . Such simulations are straightforwardly parallelized, as they feature independent trajectories. For instance, the results of figure 5 were computed from 10³ batches of 10⁵ independent simulations each; each batch takes approximately 0.8 h on 16 cores. For fixed computational resources, the observed growth of fluctuations with τ and N then leads to a trade-off between accessible time scales and system sizes. Thus, further developments will be needed for the numerical approach to be effectively applicable to late times and large systems. However, the significant suppression of fluctuations achieved using importance sampling could be particularly interesting in view of real time applications: a generalization of this method might prove useful in settings where existing techniques are significantly limited, such as higher dimensions, to explore intermediate time and system size regimes before fluctuations become dominant. The real time numerical approach [58, 59] is indeed fully analogous to the present imaginary time case, suggesting that the approach may readily generalize; however, a key difference is that for real time evolution it is not possible to make use of the τ → ∞ limit to simplify the SP equation, and a different method (e.g. recursion) will generally be needed in order to obtain the SP configuration.

In summary, in this section we used the disentanglement approach to numerically compute ground state expectation values as averages over classical stochastic trajectories. This approach is formally exact, and can be made more efficient by employing an importance sampling scheme, based on preferentially sampling trajectories close to the relevant SP configuration. Notably, this technique can be applied regardless of dimensionality. The proposed stochastic approach provides an alternative numerical method for the evaluation of ground state expectation values, which is conceptually different from existing techniques. For instance, in worldline quantum MC methods [74–77], one splits the Hamiltonian into a diagonal, 'classical' term, and a non-diagonal term; the latter is expanded in time-dependent perturbation theory, so that one is left with a sum of classical integrals which can be evaluated by MC. In contrast, here the sampling is performed over deviations from the mean field trajectory. Other classes of MC methods for spins include the previously mentioned diagrammatic approaches [69–71], based on mapping to fermions, or stochastic series expansions, whereby the whole Hamiltonian is treated perturbatively [77, 78]. The numerical performance of the stochastic approach is currently inferior to well-established numerical techniques for ground states, such as quantum MC [26, 29] or tensor network approaches [50]. However, the importance sampling method substantially extends the regime of applicability of the stochastic approach by mitigating fluctuations. This suggests that the rate of growth of fluctuations is not intrinsically fixed, and that it could be possible to suppress it even further by means of appropriate sampling schemes or approximations.

In the disentanglement approach, the Euclidean time evolution operator corresponding to the Hamiltonian (5),

$$\begin{equation}\hat{U}\left(\tau \right)\equiv {\text{e}}^{-\tau \hat{H}}={\mathrm{e}}^{\tau \left(J{\sum }_{ijab}{\mathcal{J}}_{ij}^{ab}{\hat{S}}_{i}^{a}{\hat{S}}_{j}^{b}+{\sum }_{ia}{h}_{i}^{a}{\hat{S}}_{i}^{a}\right)},\end{equation} \tag{ A1 }$$

is expressed as

$$\begin{equation}\hat{U}\left(\tau \right)=\int \mathcal{D}\varphi \enspace {\mathrm{e}}^{-{S}_{0}\left[\varphi \right]}\prod _{i}{\text{e}}^{{\xi }_{i}^{+}\left(\tau \right){\hat{S}}_{i}^{+}}\enspace {\text{e}}^{{\xi }_{i}^{z}\left(\tau \right){\hat{S}}_{i}^{z}}\enspace {\text{e}}^{{\xi }_{i}^{-}\left(\tau \right){\hat{S}}_{i}^{-}},\end{equation} \tag{ A2 }$$

where S₀ is given by equation (7) and the disentangling variables ${\xi }_{i}^{a}$ satisfy equation (9). Equation (A2) is obtained in a two-step process [24, 55, 56, 58]. First, interactions are decoupled thanks to the HS transformation [60, 61]. Following a Trotter decomposition of the exponential in equation (A1), at each time slice one has

$$\begin{equation}{\mathrm{e}}^{{\Delta}\tau J{\sum }_{ijab}{\mathcal{J}}_{ij}^{ab}{\hat{S}}_{i}^{a}{\hat{S}}_{j}^{b}}=\mathcal{C}\int \prod _{ai}\mathrm{d}{\varphi }_{i}^{a}\enspace {\mathrm{e}}^{-\frac{1}{4J}{\Delta}\tau {\sum }_{ijab}{\left({\mathcal{J}}^{-1}\right)}_{ij}^{ab}{\varphi }_{i}^{a}{\varphi }_{j}^{b}+{\Delta}\tau {\sum }_{aj}{\varphi }_{j}^{a}{\hat{S}}_{j}^{a}},\end{equation} \tag{ A3 }$$

where $\mathcal{C}$ is a normalization constant and we neglected terms O(Δτ²). Equation (A3) is an operatorial identity; the fields ${\varphi }_{i}^{a}$ are in general complex and their integration range in the complex plane is chosen in such a way as to make the integral in equation (A3) convergent [59]. It is convenient to rescale the fields as ${\varphi }_{i}^{a}\to J{\varphi }_{i}^{a}$, in order to make them dimensionless. Applying this rescaling and taking the continuum limit, equation (A3) yields

$$\begin{equation}\hat{U}\left(\tau \right)=\int \mathcal{D}\varphi \enspace {\mathrm{e}}^{-{S}_{0}\left[\varphi \right]}\mathbb{T}\prod _{i}{\mathrm{e}}^{{\sum }_{a}{\int }_{0}^{\tau }\left[{h}_{i}^{a}\left({\tau }^{\prime }\right)+J{\varphi }_{i}^{a}\left({\tau }^{\prime }\right)\right]{\hat{S}}_{i}^{a}\mathrm{d}{\tau }^{\prime }},\end{equation} \tag{ A4 }$$

where the symbol $\mathbb{T}$ denotes time ordering. Equation (A4) describes a system of non-interacting spins under the effect of complex valued stochastic fields ${\varphi }_{i}^{a}$ [24, 55]. Since interactions are decoupled inside the integral, the evolution of each spin occurs over its (complexified) Bloch sphere. Time-ordered exponentials can then be expressed in terms of ordinary exponentials by means of a Lie-algebraic disentanglement transformation, also known as Wei–Norman–Kolokolov transformation [24, 55, 56, 62, 63]. Namely, at each lattice site one has

$$\begin{equation}\mathbb{T}\enspace {\mathrm{e}}^{{\sum }_{a}{\int }_{0}^{\tau }\left[{h}_{i}^{a}\left({\tau }^{\prime }\right)+J{\varphi }_{i}^{a}\left({\tau }^{\prime }\right)\right]{\hat{S}}_{i}^{a}\mathrm{d}{\tau }^{\prime }}={\text{e}}^{{\xi }_{i}^{+}\left(\tau \right){\hat{S}}_{i}^{+}}\enspace {\text{e}}^{{\xi }_{i}^{z}\left(\tau \right){\hat{S}}_{i}^{z}}\enspace {\text{e}}^{{\xi }_{i}^{-}\left(\tau \right){\hat{S}}_{i}^{-}}.\end{equation} \tag{ A5 }$$

This amounts to parameterizing the trajectory of each spin on its Bloch sphere in terms of a set of coordinates ${\xi }_{i}^{a}$, termed the disentangling variables. Equation (A5) can be seen as the defining equation of ${\xi }_{i}^{a}$; differentiating both sides of (A5) and equating the coefficients that multiply the spin operators yields the SDEs (9). Alternatively, these can be obtained from differential geometry [24]. The initial conditions ${\xi }_{i}^{a}\left(0\right)=0$ are fixed by the requirement $\hat{U}\left(0\right)=\mathbb{1}$. The discussion of this section can be readily generalized to the modified time evolution operator $\hat{\mathcal{U}}\left(\tau \right)$, given by equation (2); the initial conditions of the disentangling variables are then given by (10).

While the formalism outlined in the previous section is fully general, in this section we show in greater detail how the disentanglement transformation works in higher dimensional settings of particular physical interest, providing useful formulae for analytical and numerical applications. Let us consider a Hamiltonian describing a system on a D-dimensional hypercubic lattice:

$$\begin{equation}\hat{H}=-\sum _{\boldsymbol{i}\boldsymbol{j}}\sum _{ab}{\mathcal{J}}_{\boldsymbol{i}\boldsymbol{j}}^{ab}{\hat{S}}_{\boldsymbol{i}}^{a}{\hat{S}}_{\boldsymbol{j}}^{b}-\sum _{\boldsymbol{i}}\sum _{a}{h}_{\boldsymbol{i}}^{a}{\hat{S}}_{\boldsymbol{i}}^{a}.\end{equation} \tag{ A6 }$$

We focus on the case where $\mathcal{J}$ only couples spins along the lattice axes, i.e. the sites i , j coupled by ${\mathcal{J}}_{\boldsymbol{i}\boldsymbol{j}}$ differ by a single index, i_d ≠ j_d . For this choice of $\mathcal{J}$, one has

$$\begin{equation}{\mathcal{J}}_{\boldsymbol{i}\boldsymbol{j}}^{ab}=\sum _{d}{J}_{d}{\left({\mathcal{J}}^{d}\right)}_{{i}_{d}{j}_{d}}^{ab}\prod _{{d}^{\prime }\ne d}{\delta }_{{i}_{{d}^{\prime }}{j}_{{d}^{\prime }}},\end{equation} \tag{ A7 }$$

where J_d are interaction strengths and the matrices ${\mathcal{J}}^{d}$ couple spins along the dimension d. ${\mathcal{J}}^{d}$ can be seen as 3N_d × 3N_d matrices ${\mathcal{J}}_{\alpha \beta }^{d}$ by introducing multi-component indices α = {i_d , a}, β = {j_d , b}. One has

$$\begin{equation}\sum _{\boldsymbol{i}\boldsymbol{j}}\sum _{ab}{\mathcal{J}}_{\boldsymbol{i}\boldsymbol{j}}^{ab}{\hat{S}}_{\boldsymbol{i}}^{a}{\hat{S}}_{\boldsymbol{j}}^{b}=\sum _{d=1}^{D}{J}_{d}\sum _{\boldsymbol{i}}{\hat{S}}_{\boldsymbol{i}}^{a}\sum _{{j}_{d}}{\left({\mathcal{J}}^{d}\right)}_{{i}_{d}{j}_{d}}^{ab}{\hat{S}}_{{i}_{1}\dots {j}_{d}\dots {i}_{D}}^{b}.\end{equation} \tag{ A8 }$$

Exponentiating the interaction term in (A6) and considering an infinitesimal time slice, one obtains

$$\begin{equation}{\mathrm{e}}^{{\Delta}\tau {\sum }_{\boldsymbol{i}\boldsymbol{j}ab}{\mathcal{J}}_{\boldsymbol{i}\boldsymbol{j}}^{ab}{\hat{S}}_{\boldsymbol{i}}^{a}{\hat{S}}_{\boldsymbol{j}}^{b}}=\prod _{d=1}^{D}\prod _{\boldsymbol{i}\ne {i}_{d}}{\mathrm{e}}^{{\Delta}\tau {\sum }_{{i}_{d}{j}_{d}ab}{J}_{d}{\left({\mathcal{J}}^{d}\right)}_{{i}_{d}{j}_{d}}^{ab}{\hat{S}}_{\boldsymbol{i}}^{a}{\hat{S}}_{{i}_{1}\dots {j}_{d}\dots {i}_{D}}^{b}}.\end{equation} \tag{ A9 }$$

We can apply the HS transformation to each term:

$$\begin{align}\hfill & {\mathrm{e}}^{{\Delta}\tau {\sum }_{{i}_{d}{j}_{d}ab}{J}_{d}{\left({\mathcal{J}}^{d}\right)}_{{i}_{d}{j}_{d}}^{ab}{\hat{S}}_{\boldsymbol{i}}^{a}{\hat{S}}_{{i}_{1}\dots {j}_{d}\dots {i}_{D}}^{b}}=\mathcal{C}\int \prod _{{i}_{d}=1}^{{N}_{d}}\left(\mathrm{d}{\varphi }_{\boldsymbol{i}}^{d}\right)\hfill \\ \hfill & \quad {\times}{\mathrm{e}}^{-\frac{{\Delta}\tau }{4{J}_{d}}{\sum }_{{i}_{d}{j}_{d}ab}{\left[{\left({\mathcal{J}}^{d}\right)}^{-1}\right]}_{{i}_{d}{j}_{d}}^{ab}{\left({\varphi }^{d}\right)}_{\boldsymbol{i}}^{a}{\left({\varphi }^{d}\right)}_{{i}_{1}\dots {j}_{d}\dots {i}_{D}}^{b}+{\Delta}\tau {\sum }_{{j}_{d}a}{\left({\varphi }^{d}\right)}_{{i}_{1}\dots {j}_{d}\dots {i}_{D}}^{a}{\hat{S}}_{{i}_{1}\dots {j}_{d}\dots {i}_{D}}^{a}}.\hfill \end{align} \tag{ A10 }$$

Taking the continuum limit and rescaling φ^d → J_d φ^d , this yields

$$\begin{equation}{\text{e}}^{-\tau \hat{H}}=\int \mathcal{D}\varphi \enspace {\mathrm{e}}^{-{S}_{0}\left[\varphi \right]}\mathbb{T}\enspace {\mathrm{e}}^{{\int }_{0}^{\tau }\mathrm{d}\tau \left[\sum _{d}{J}_{d}\sum _{\boldsymbol{i}a}{\left({\varphi }^{d}\right)}_{\boldsymbol{i}}^{a}+{h}_{\boldsymbol{i}}^{a}\right]{\hat{S}}_{\boldsymbol{i}}^{a}}\end{equation} \tag{ A11 }$$

where the noise action in D dimensions is given by

$$\begin{equation}{S}_{0}\left[\varphi \right]=\frac{1}{4}{\int }_{0}^{\tau }\mathrm{d}\tau \sum _{d}{J}_{d}\sum _{\boldsymbol{i}{j}_{d}ab}{\left[{\left({\mathcal{J}}^{d}\right)}^{-1}\right]}_{{i}_{d}{j}_{d}}^{ab}{\left({\varphi }^{d}\right)}_{\boldsymbol{i}}^{a}{\left({\varphi }^{d}\right)}_{{i}_{1}\dots {j}_{d}{i}_{D}}^{b}.\end{equation} \tag{ A12 }$$

It can be seen that for a D dimensional system with N = N₁ ×⋯× N_d spins, one needs in general to introduce 3DN HS fields. The individual time ordered exponentials in equation (A11) can then be expressed in terms of ordinary exponentials, as done in appendix A.1. As in the 1D case, a change of variables can be performed to make the noise action S₀ diagonal. Let us introduce the notation ${\bar{\boldsymbol{i}}}_{d}\equiv {i}_{1}\dots {i}_{d-1}{i}_{d+1}\dots {i}_{D}$ and ${\left({\varphi }^{d}\right)}_{\boldsymbol{i}}^{a}={\left({\varphi }^{d}\right)}_{{\bar{\boldsymbol{i}}}_{d}\alpha }$ with α = {i_d , a}. Equation (A12) is then diagonalized by the transformation

$$\begin{equation}{\left({\varphi }^{d}\right)}_{\boldsymbol{i}}^{a}={\left({\varphi }^{d}\right)}_{{\bar{\boldsymbol{i}}}_{d}\alpha }=\sum _{\beta =1}^{3{N}_{d}}{\left({O}^{d}\right)}_{\alpha \beta }{\left({\phi }^{d}\right)}_{{\bar{\boldsymbol{i}}}_{d}\beta },\end{equation} \tag{ A13 }$$

where O^d is a 3N_d × 3N_d matrix defined as for the 1D case, but in terms of ${\mathcal{J}}^{d}$. Using equation (A13), we obtain

$$\begin{equation}{\text{e}}^{\tau {\sum }_{\boldsymbol{i}\boldsymbol{j}}{\mathcal{J}}_{\boldsymbol{i}\boldsymbol{j}}{\hat{S}}_{\boldsymbol{i}}^{z}{\hat{S}}_{\boldsymbol{j}}^{z}}=\int \mathcal{D}\phi \enspace {\mathrm{e}}^{-{S}_{0}\left[\phi \right]+{\int }_{0}^{\tau }\mathrm{d}\tau {\sum }_{\boldsymbol{i}a}{\hat{S}}_{\boldsymbol{i}}^{a}{\sum }_{d}{J}_{d}\sum _{\beta =1}^{3{N}_{d}}{\left({O}^{d}\right)}_{\alpha \beta }{\left({\phi }^{d}\right)}_{{\bar{\boldsymbol{i}}}_{d}\beta }}\end{equation} \tag{ A14 }$$

with

$$\begin{equation}{S}_{0}\left[\phi \right]=\frac{1}{2}{\int }_{0}^{{\tau }_{f}}\mathrm{d}\tau \sum _{d}\sum _{\boldsymbol{i}a}{\left[{\left({\phi }^{d}\right)}_{\boldsymbol{i}}^{a}\right]}^{2}.\end{equation} \tag{ A15 }$$

The stochastic representation (12) of the Euclidean time evolution operator is formally exact; this implies that the statistics of the classical disentangling variables $\xi =\left\{{\xi }_{i}^{a}\right\}$ contain all the information about the corresponding quantum problem. In the case of real time evolution, this observation was drawn upon in [58, 59] to numerically investigate the relation between fluctuations in the disentangling variables and dynamical quantum phase transitions (DQPTs) [85, 86]. Here we consider the imaginary time behavior of the disentangling variables, which encodes all information about the ground state of the corresponding quantum problem. For definiteness, we consider the quantum Ising model, given by the Hamiltonian (20). For the one-dimensional quantum Ising chain, the general result (9) specializes to the Euclidean Ising SDEs [24, 58, 59]

$$\begin{equation}{\dot {\xi }}_{i}^{+}\left(\tau \right)=\frac{{\Gamma}}{2}\left(1-{{\xi }_{i}^{+}}^{2}\right)+J{\xi }_{i}^{+}\sum _{j}{O}_{ij}{\phi }_{j},\end{equation} \tag{ B1a }$$

$$\begin{equation}{\dot {\xi }}_{i}^{z}\left(\tau \right)=-{\Gamma}{\xi }_{i}^{+}+J\sum _{j}{O}_{ij}{\phi }_{j},\end{equation} \tag{ B1b }$$

$$\begin{equation}{\dot {\xi }}_{i}^{-}\left(\tau \right)=\frac{{\Gamma}}{2}\enspace \mathrm{exp}\enspace {\xi }_{i}^{z},\end{equation} \tag{ B1c }$$

which are here expressed in terms of the fields ϕ_i that diagonalize the noise action (7).

To the best of our current knowledge, equation (B1) are only exactly solvable in the classical (Γ = 0) and non-interacting (J = 0) cases. This was discussed in [59] for real time evolution and in the special case ${\xi }_{i}^{a}\left(0\right)=0$; here we consider Euclidean time and general initial conditions, as it is relevant for our current purposes. For the present discussion, we focus on the one-dimensional case, which is sufficient to illustrate the relevant properties of the disentangling variables; the higher-dimensional version of equation (B1) is given by equation (C8) in appendix C. In the classical case with Γ = 0, the non-linear term in equation (B1a) vanishes and ${\xi }_{i}^{+}$ performs driftless geometric Brownian motion. This is exactly solvable, giving

$$\begin{equation}{\xi }_{i}^{+}\left(\tau \right)={\xi }_{i}^{+}\left(0\right)\enspace \mathrm{exp}\left[\sum _{j}{O}_{ij}{\int }_{0}^{\tau }{\phi }_{j}\left(s\right)\mathrm{d}s\right],\end{equation} \tag{ B2 }$$

where we used ${\left(O{O}^{\mathrm{T}}\right)}_{ii}\propto {\mathcal{J}}_{ii}=0$. In the classical limit, ${\xi }_{i}^{z}$ is decoupled from ${\xi }_{i}^{+}$ and satisfies Brownian motion:

$$\begin{equation}{\xi }_{i}^{z}\left(\tau \right)={\xi }_{i}^{z}\left(0\right)+{\int }_{0}^{\tau }\sum _{j}{O}_{ij}{\phi }_{j}\left(s\right)\mathrm{d}s,\end{equation} \tag{ B3 }$$

while ${\xi }_{i}^{-}\left(\tau \right)={\xi }_{i}^{-}\left(0\right)$.

In the non-interacting limit J = 0, equation (B1) become deterministic and solvable, yielding

$$\begin{equation}{\xi }_{i}^{+}\left(\tau \right)={\xi }_{i}^{+}\left(0\right)+\frac{1-{\xi }_{i}^{+2}\left(0\right)}{{\xi }_{i}^{+}\left(0\right)+\mathrm{coth}\left({\Gamma}\tau /2\right)},\end{equation} \tag{ B4a }$$

$$\begin{equation}{\xi }_{i}^{z}\left(\tau \right)={\xi }_{i}^{z}\left(0\right)-2\enspace \mathrm{log}\left[\mathrm{cosh}\left({\Gamma}\tau /2\right)+{\xi }_{i}^{+}\left(0\right)\mathrm{sinh}\left({\Gamma}\tau /2\right)\right],\end{equation} \tag{ B4b }$$

$$\begin{equation}{\xi }_{i}^{-}\left(\tau \right)={\xi }_{i}^{-}\left(0\right)+\frac{\mathrm{exp}\left[{\xi }_{i}^{z}\left(0\right)\right]}{{\xi }_{i}^{+}\left(0\right)+\mathrm{coth}\left({\Gamma}\tau /2\right)}.\end{equation} \tag{ B4c }$$

In the general case with finite Γ and J, the SDEs (B1) cannot be solved exactly to the best of our knowledge. However, analytical insights about (B1) can still be obtained. Of particular interest is the behavior of the variables ${\xi }_{i}^{+}$: as observed in [24, 58, 59], these play a key role, being the primary source of non-linearity in (B1) (the variable ${\xi }_{i}^{-}$ is seldom needed to compute observables) and the sole disentangling variable whose equation of motion is autonomous, not involving any other ${\xi }_{i}^{a}$. The stationary probability distribution attained at late times by ${\xi }_{i}^{+}\left(\tau \right)$ was obtained in [24, 64]. Additional information is encoded in the moment-generating function G_i (λ, τ) of ${\xi }_{i}^{+}\left(\tau \right)$, which satisfies ${\partial }_{\lambda }^{n}{G}_{i}\left(\lambda ,\tau \right){\vert }_{\lambda =0}={\langle {\xi }_{i}^{+n}\left(\tau \right)\rangle }_{\phi }$ and gives access to the Euclidean time-dependent moments of ${\xi }_{i}^{+}$. To compute this, we define the stochastic function ${g}_{i}\left(\lambda ,\tau \right)\equiv {\mathrm{e}}^{\lambda {\xi }_{i}^{+}\left(\tau \right)}$, such that ${G}_{i}\left(\lambda ,\tau \right)\equiv {\langle {g}_{i}\left(\lambda ,\tau \right)\rangle }_{\phi }$. The equation of motion of g_i (τ) is obtained by applying the Ito chain rule [57, 87]:

$$\begin{align}\hfill \frac{\mathrm{d}}{\mathrm{d}\tau }{g}_{i}\left(\lambda ,\tau \right)& =\lambda \frac{{\Gamma}}{2}\left(1-{{\xi }_{i}^{+}}^{2}\right){g}_{i}\left(\lambda ,\tau \right)+\lambda {\xi }_{i}^{+}\sum _{j}{O}_{ij}{\phi }_{j}{g}_{i}\left(\lambda ,\tau \right)\hfill \\ \hfill & \quad +\frac{1}{2}{\lambda }^{2}{{\xi }_{i}^{+}}^{2}\sum _{j}{O}_{ij}{O}_{ij}{g}_{i}\left(\lambda ,\tau \right).\hfill \end{align} \tag{ B5 }$$

It can be easily shown that the matrix OO^T is proportional to $\mathcal{J}$ and hence has no diagonal term [59]; this implies that the Ito drift term proportional to ∑_j O_ij O_ij gives no contribution. It is also convenient to write ${\xi }_{i}^{+}{g}_{i}=\frac{\partial }{\partial \lambda }{g}_{i}$. With these simplifications, we obtain

$$\begin{equation}\frac{\mathrm{d}}{\mathrm{d}\tau }{g}_{i}\left(\lambda ,\tau \right)=\left[\lambda \frac{{\Gamma}}{2}\left(1-\frac{{\partial }^{2}}{\partial {\lambda }^{2}}\right)+\sum _{j}{O}_{ij}{\phi }_{j}\frac{\partial }{\partial \lambda }\right]{g}_{i}\left(\lambda ,\tau \right).\end{equation} \tag{ B6 }$$

Considering the expectation value of equation (B6) and using the property of Ito calculus ${\langle {g}_{i}\left(\tau \right){\phi }_{j}\left(\tau \right)\rangle }_{\phi }=0\enspace \enspace \forall \enspace i,j$ we obtain the partial differential equation satisfied by the moment-generating function:

$$\begin{equation}\frac{\partial }{\partial t}{G}_{i}\left(\lambda ,\tau \right)=\left[\lambda \frac{{\Gamma}}{2}\left(1-\frac{{\partial }^{2}}{\partial {\lambda }^{2}}\right)\right]{G}_{i}\left(\lambda ,\tau \right),\end{equation} \tag{ B7 }$$

with initial conditions G_i (0, τ) = G_i (λ, 0) = 1. Equation (B7) can be solved exactly, yielding

$$\begin{equation}{G}_{i}\left(\lambda ,\tau \right)=\mathrm{exp}\left[\lambda \left({\xi }_{i}^{+}\left(0\right)+\frac{1-{\xi }_{i}^{+2}\left(0\right)}{{\xi }_{i}^{+}\left(0\right)+\mathrm{coth}\left({\Gamma}\tau /2\right)}\right)\right].\end{equation} \tag{ B8 }$$

This result predicts that all moments of ${\xi }_{i}^{+}\left(\tau \right)$ are given by powers of the deterministic trajectory obtained in the non-interacting case with J = 0: the moments of each individual ${\xi }_{i}^{+}$ contain no information about the interacting quantum system. All information is therefore encoded in the correlations between variables at different sites. We note that the findings of the present section do not apply to real time evolution: in that case, the moments of ${\xi }_{i}^{+}\left(t\right)$ for non-zero J differ from the non-interacting result. This discrepancy can be traced back to the failure of the analytic continuation of equation (B8) to real time.

Since the information about interactions is contained in the joint statistics of the ${\xi }_{i}^{+}$ variables, we investigate the joint probability distribution $P\left[{\xi }^{+}\right]\equiv P\left[\left\{{\xi }_{i}^{+}\right\}\right]$. The stochastic process ${\xi }_{i}^{+}$ has drift and diffusion

$$\begin{equation}{a}_{i}\left({\xi }_{i}^{+}\right)=\frac{{\Gamma}}{2}\left(1-{\xi }_{i}^{+2}\right),\end{equation} \tag{ B9a }$$

$$\begin{equation}{B}_{ij}\left({\xi }_{i}^{+}\right)={\xi }_{i}^{+}{O}_{ij}\end{equation} \tag{ B9b }$$

respectively. The probability distribution of its realizations is given by [24, 88–90]

$$\begin{equation}P\left[{\xi }^{+}\right]={\mathcal{C}}_{\xi }\enspace {\mathrm{e}}^{-I\left[{\xi }^{+}\right]},\end{equation} \tag{ B10 }$$

where ${\mathcal{C}}_{\xi }$ is a normalization constant and

$$\begin{equation}I\left[{\xi }^{+}\right]={\int }_{0}^{\tau }\mathrm{d}{\tau }^{\prime }L\left({\xi }^{+},{\dot {\xi }}^{+}\right),\end{equation} \tag{ B11 }$$

$$\begin{equation}L\left({\xi }^{+},{\dot {\xi }}^{+}\right)=\frac{1}{2}\sum _{ij}\left[{\dot {\xi }}_{i}^{+}-{a}_{i}\left({\xi }_{i}^{+}\right)\right]{\mathcal{B}}_{ij}^{-1}\left({\xi }^{+}\right)\left[{\dot {\xi }}_{j}-{a}_{j}\left({\xi }_{j}^{+}\right)\right],\end{equation} \tag{ B12 }$$

$$\begin{equation}{\mathcal{B}}_{ij}\left({\xi }^{+}\right)=\sum _{k}{B}_{ik}\left({\xi }_{i}^{+}\right){B}_{jk}\left({\xi }_{j}^{+}\right).\end{equation} \tag{ B13 }$$

Equation (B9) give ${\mathcal{B}}_{ij}\left({\xi }^{+}\right)=2J{\mathcal{J}}_{ij}{\xi }_{i}^{+}{\xi }_{j}^{+}$ and

$$\begin{equation}L\left({\xi }^{+},{\dot {\xi }}^{+}\right)=\frac{1}{4J}\sum _{ij}\frac{1}{{\xi }_{i}^{+}{\xi }_{j}^{+}}\left[{\dot {\xi }}_{i}^{+}-\frac{{\Gamma}}{2}\left(1-{{\xi }_{i}^{+}}^{2}\right)\right]{\mathcal{J}}_{ij}^{-1}\left[{\dot {\xi }}_{j}^{+}-\frac{{\Gamma}}{2}\left(1-{{\xi }_{j}^{+}}^{2}\right)\right].\end{equation} \tag{ B14 }$$

Equation (B10) provides the measure when the stochastic expression for an observable is expressed as a path integral over the variables ${\xi }_{i}^{+}$ rather than the fields ϕ_i [24]. When sampling according to the distribution (B10), the likeliest trajectory is obtained by extremizing the weight I[ξ⁺] with respect to ${\xi }_{i}^{+}\left(\tau \right)$. By solving the corresponding Euler–Lagrange equations, we readily see that the dominant trajectory is the non-interacting solution (B4a). Therefore, when applying the stochastic approach using direct sampling [58, 59], one typically samples trajectories which are nearly non-interacting.

We note that in the large Γ limit equation (B10) takes a large deviation form [91]. Since Γ multiplies time in ${\xi }_{\text{NI}}^{+}$, we rescale time as $\tilde {\tau }=\tau {\Gamma}$. The corresponding rescaled stochastic equation for ${\xi }_{i}^{+}$ is

$$\begin{equation}{\dot {\xi }}_{i}^{+}\left(\tilde {\tau }\right)=\frac{1}{2}\left(1-{{\xi }_{i}^{+}}^{2}\right)+{\epsilon}{\xi }_{i}^{+}\sum _{j}{O}_{ij}{\phi }_{j},\end{equation} \tag{ B15 }$$

where we have defined the noise strength ≡ 1/Γ. The limit Γ → ∞ is therefore equivalent to the small-noise limit of (B15). SDEs in the limit of small noise are described by the Freidlin–Wentzell large deviation theory [92, 93]: ${\xi }_{i}^{+}$ obeys a LDP, with rate ⁻² and rate function $\mathcal{I}\left[{\xi }^{+}\right]\equiv {{\epsilon}}^{2}I\left[{\xi }^{+}\right]$. In this small- limit, the trajectories ${\xi }_{i}^{+}$ are approximately Gaussian distributed around the likeliest trajectory ${\xi }_{\text{NI}}^{+}$ [91]:

$$\begin{equation}P\left[{\xi }_{i}^{+}\right]\sim {\mathrm{e}}^{-\frac{{{\epsilon}}^{-2}}{2}{\int }_{0}^{\tilde {\tau }}{\sum }_{ij}{\mathcal{I}}_{ij}^{\left(2\right)}\left[{\xi }_{i}^{+}\left({\tilde {\tau }}^{\prime }\right)-{\xi }_{\mathrm{N}\mathrm{I}}^{+}\left({\tilde {\tau }}^{\prime }\right)\right]\left[{\xi }_{j}^{+}\left({\tilde {\tau }}^{\prime }\right)-{\xi }_{\mathrm{N}\mathrm{I}}^{+}\left({\tilde {\tau }}^{\prime }\right)\right]\mathrm{d}{\tilde {\tau }}^{\prime }},\end{equation} \tag{ B16 }$$

where the second variation ${\mathcal{I}}_{ij}^{\left(2\right)}$ is given by

$$\begin{equation}{\mathcal{I}}_{ij}^{\left(2\right)}\equiv {\left.\frac{{\delta }^{2}\mathcal{I}}{\delta {\xi }_{i}^{+}\left({\tilde {\tau }}^{\prime }\right)\delta {\xi }_{j}^{+}\left({\tilde {\tau }}^{\prime }\right)}\right\vert }_{{\xi }_{\mathrm{N}\mathrm{I}}^{+}}.\end{equation} \tag{ B17 }$$

Thus, trajectories that deviate significantly from the non-interacting limit are exponentially suppressed. The large deviation formalism also applies to real time evolution, where again the dominant trajectory is given by the deterministic result ${\xi }_{\text{NI}}^{+}\left(t\right)$. However, in contrast to ${\xi }_{\text{NI}}^{+}\left(\tau \right)$, ${\xi }_{\text{NI}}^{+}\left(t\right)$ has an infinite number of singularities as a function of time [59]. This leads to a breakdown of the expansion about the non-interacting SP, which can be expected to have consequences for sampling. Even for large Γ, regions in time that are close to the singularities in the SP trajectory are expected to be associated with enhanced fluctuations, leading to difficulties in sampling. This observation may lie at the root of the enhanced fluctuations of the disentangling variables found in the vicinity of DQPTs [85, 86], reported in [58, 59].

The SP equation (30) for the Loschmidt amplitude can be solved recursively, exploiting the intuition that the SP field configuration φ_SP(τ') ≡ φ_SP(τ'|τ_f ) should change little if τ_f is increased by a small amount Δt. In practice, one assumes

$$\begin{equation}{\varphi }_{\text{SP}}\left({\tau }^{\prime }\vert {\tau }_{f}+{\Delta}t\right)\approx {\varphi }_{\text{SP}}\left({\tau }^{\prime }\vert {\tau }_{f}\right)\end{equation} \tag{ C1 }$$

for τ' < τ_f + Δt. The field φ_SP(τ'|τ_f ) is then used to compute ${\xi }_{i}^{+}{\vert }_{\text{SP}}$ and Ξ_ij |_SP. Using these quantities, one can in turn produce a better approximation of φ_SP(τ'|τ_f + Δt) according to the SP equation (30). This procedure can be iterated until the field configuration has converged to a desired level of accuracy. The convergence of the recursion is determined by defining a quantity ɛ which measures how much the approximate SP field varies after an iteration of the algorithm. A suitable definition is

$$\begin{equation}\varepsilon \equiv \frac{1}{k}\sum _{m=1}^{k}\left\vert {\bar{\varphi }}_{\text{SP}}\left({\tau }_{m}\vert {\tau }_{f}+{\Delta}\tau \right)-{\varphi }_{\text{SP}}\left({\tau }_{m}\vert {\tau }_{f}+{\Delta}\tau \right)\right\vert \end{equation} \tag{ C2 }$$

where φ_SP and ${\bar{\varphi }}_{\text{SP}}$ are the old and updated estimates of the SP field respectively, evaluated at the discrete times τ_m . Convergence is then defined as ɛ < ɛ*, where ɛ* is a threshold of choice. The runtime of this recursive algorithm scales quadratically with the number of time steps n; this is because for each 1 < k < n one needs to perform k calculations in order to compute ${\xi }_{i}^{+}{\vert }_{\text{SP}}$, so that summing over all k the total number of calculations to perform is of order n(n + 1)/2. In principle, the computational cost is further increased by having to repeat each step multiple times to attain convergence. However, for reasonable values of the threshold ɛ*, numerical evaluation shows that the recursive algorithm has rapid convergence, typically requiring only 1–2 iterations. From recursively solving the SP equation, we find that for sufficiently large τ_f the value φ_SP(τ'|τ_f ) with τ' ≪ τ_f no longer changes with τ_f , settling to a value φ_SP(τ'|∞) ≡ φ_SP(τ'); this is illustrated in figure 10(a). Because of this, when recursively solving the SP equation one only needs to update the SP configuration at the times τ' such that φ_SP(τ'|τ_f ) ≠ φ_SP(τ') to a desired level of precision; this speeds up the recursive solution significantly. The SP equation (30) prescribes that the value of the SP field at the end time is always φ_SP(τ_f |τ_f ) = −1. Thus, the SP field φ_SP(τ'|τ_f ) cannot attain a steady state, i.e. for finite τ_f there exists no time scale τ_SS such that ∂_τ' φ_SP(τ', τ_f ) ≈ 0 ∀ τ' > τ_SS. However, the numerical results show that for sufficiently large τ_f the SP field φ_SP(τ', τ_f ) attains a plateau value at times 0 ≪ τ' ≪ τ_f ; this is illustrated in figure 10(b). The extent of the plateau grows as τ_f is increased; since the action is extensive in time, the plateau value provides the dominant contribution to observables in the large τ_f limit. The plateau value of φ_SP can be found analytically, as discussed in section 2.3; the analytical results are in perfect agreement with the numerical solution.

Zoom In Zoom Out Reset image size

In the disentanglement formalism, the Euclidean time evolution of an observable $\mathcal{O}$ is given by

$$\begin{equation}\langle {\psi }_{0}\vert \hat{U}\left(\tau \right)\hat{\mathcal{O}}\hat{U}\left(\tau \right)\vert {\psi }_{0}\rangle =\int \mathcal{D}\phi \enspace {\mathrm{e}}^{-{S}_{0}\left[\phi \right]}{F}_{\mathcal{O}}\left[\phi \right],\end{equation} \tag{ C3 }$$

where $\phi =\left\{{\phi }_{\mathrm{f},i}^{a},{\phi }_{\mathrm{b},i}^{a}\right\}$ collectively denotes the two sets of HS fields introduced to decouple the two time-evolution operators, and the classical function ${F}_{\mathcal{O}}$ is given by equation (14). As discussed in section 2.2, the trajectory yielding the largest contribution to the integral can be found by extremizing the effective action

$$\begin{equation}{S}_{\mathcal{O}}\equiv {S}_{0}\left[\phi \right]-\mathrm{log}\enspace {f}_{\mathcal{O}}\left[\phi \right].\end{equation} \tag{ C4 }$$

Consider the normalization function, corresponding to setting $\hat{\mathcal{O}}=\mathbb{1}$ in (C3) and given by equation (15). For the 1D quantum Ising model, the effective action for this quantity is given by

$$\begin{align}\hfill {S}_{\mathbb{1}}& =\frac{1}{2}{\int }_{0}^{{\tau }_{f}}\mathrm{d}\tau \sum _{i}\left[\frac{J}{2}\sum _{j}{\mathcal{J}}_{ij}^{-1}\left[{\varphi }_{\mathrm{f},i}{\varphi }_{\mathrm{f},j}+{\varphi }_{\mathrm{b},i}^{{\ast}}{\varphi }_{\mathrm{b},j}^{{\ast}}\right]-{\Gamma}{\xi }_{\mathrm{f},i}^{+}-{\Gamma}{\xi }_{\mathrm{b},i}^{+{\ast}}+J{\varphi }_{\mathrm{f},i}+J{\varphi }_{\mathrm{b},i}^{{\ast}}\right]\hfill \\ \hfill & \quad -\sum _{i}\mathrm{log}\left[1+{\xi }_{\mathrm{f},i}^{+}\left({\tau }_{f}\right){\xi }_{\mathrm{b},i}^{+{\ast}}\left({\tau }_{f}\right)\right].\hfill \end{align} \tag{ C5 }$$

By varying equation (C5), we obtain the SP equations for the normalization:

$$\begin{equation}J\sum _{j}{\mathcal{J}}_{ij}^{-1}{\varphi }_{\mathrm{f},j}\left({\tau }^{\prime }\right){\vert }_{\text{SP}}={\Gamma}{\int }_{0}^{{\tau }_{f}}{{\Xi}}_{\mathrm{f},i}\left(\tau ,{\tau }^{\prime }\right)\mathrm{d}\tau {\vert }_{\text{SP}}-J+{\left.\frac{2{\xi }_{\mathrm{b},i}^{+{\ast}}{{\Xi}}_{i}\left({\tau }_{f},{\tau }^{\prime }\right)}{1+{\xi }_{\mathrm{f},i}^{+}\left({\tau }_{f}\right){\xi }_{\mathrm{b},i}^{+{\ast}}\left({\tau }_{f}\right)}\right\vert }_{\text{SP}}.\end{equation} \tag{ C6 }$$

The same equation is satisfied by φ_b,i |_SP, with the replacement f ↔ b. By direct substitution, one readily verifies that the plateau of the Loschmidt amplitude SP is a fixed point of equation (C6) at all times. Thus, choosing the MF ground state as the initial state eliminates both the transient and the late-time behavior of equation (C6) (in contrast, for any initial state, the solution of equation (30) deviates from the plateau at late times due to the boundary condition ϕ(τ_f ) = −1). More generally, local observables expressed in a translationally invariant way correspond to stochastic functions

$$\begin{equation}{f}_{\mathcal{O}}\left({\tau }_{f}\right)={f}_{\mathbb{1}}\left({\tau }_{f}\right)\sum _{i}{\bar{f}}_{\mathcal{O},i}\left({\tau }_{f}\right),\end{equation} \tag{ C7 }$$

where ${\bar{f}}_{\mathcal{O},i}\left({\tau }_{f}\right)$ is a function of ${\xi }_{\mathrm{f},i}^{+}\left({\tau }_{f}\right)$ only. For instance, for the magnetization one has ${\bar{f}}_{\mathcal{M},i}=\left(1-{\xi }_{\mathrm{f},i}^{+}{\xi }_{\mathrm{b},i}^{+{\ast}}\right)/\left(1+{\xi }_{\mathrm{f},i}^{+}{\xi }_{\mathrm{b},i}^{+{\ast}}\right)$; see equation (16). It can be readily seen that the SP equation obtained by extremizing the effective action for (C7) differs by equation (C6) by a term proportional to 1/N. Furthermore, the extra term is also proportional to Ξ_SP(τ_f , τ); at the plateau, one has ${{\Xi}}_{\mathrm{P}}\left({\tau }_{f},\tau \right)\propto {\mathrm{e}}^{-\left({\Gamma}{\xi }_{\mathrm{P}}^{+}-J{\phi }_{\mathrm{P}}\right)\left({\tau }_{f}-\tau \right)}$, so that the extra term is inconsequential as τ_f → ∞. Thus, the SP equation for any local observable differs from (C6) by a term which is suppressed both as τ_f → ∞ and as N → ∞. This implies that the plateau SP trajectory for the Loschmidt amplitude can be used for all other ground state expectation values, both for analytical and numerical applications.

For the D-dimensional quantum Ising model, the Euclidean SDEs are given by

$$\begin{equation}{\dot {\xi }}_{\boldsymbol{i}}^{+}=\frac{{\Gamma}}{2}\left(1-{{\xi }_{\boldsymbol{i}}^{+}}^{2}\right)+{\xi }_{\boldsymbol{i}}^{+}\sum _{d=1}^{D}{J}_{d}{\varphi }_{\boldsymbol{i}}^{d},\end{equation} \tag{ C8a }$$

$$\begin{equation}{\dot {\xi }}_{\boldsymbol{i}}^{z}=-{\Gamma}{\xi }_{\boldsymbol{i}}^{+}+\sum _{d=1}^{D}{J}_{d}{\varphi }_{\boldsymbol{i}}^{d},\end{equation} \tag{ C8b }$$

$$\begin{equation}{\dot {\xi }}_{\boldsymbol{i}}^{-}=\frac{{\Gamma}}{2}\enspace \mathrm{exp}\enspace {\xi }_{\boldsymbol{i}}^{z},\end{equation} \tag{ C8c }$$

with multicomponent indices i = {i₁, ..., i_D }. The Euclidean Loschmidt amplitude is given by

$$\begin{equation}A\left({\tau }_{f}\right)=\int \mathcal{D}\varphi \enspace {\mathrm{e}}^{-{S}_{0}\left[\varphi \right]-\frac{1}{2}{\int }_{0}^{{\tau }_{f}}{\sum }_{\boldsymbol{i}}\mathrm{d}\tau \left[\sum _{d=1}^{D}{J}_{d}{\varphi }_{\boldsymbol{i}}^{d}\left(\tau \right)-{\Gamma}{\xi }_{\boldsymbol{i}}^{+}\left(\tau \right)\right]}.\end{equation} \tag{ C9 }$$

The SP equation obtained by varying the effective action with respect to ${\varphi }_{\boldsymbol{i}}^{d}\left({\tau }^{\prime }\right)$ is then given by

$$\begin{equation}{\varphi }_{\boldsymbol{i}}^{d}\left({\tau }^{\prime }\right){\vert }_{\text{SP}}=\frac{{\Gamma}}{{J}_{d}}\sum _{{j}_{d}}{\mathcal{J}}_{{i}_{d}{j}_{d}}^{d}{\int }_{0}^{{\tau }_{f}}{\left.{{\Xi}}_{{i}_{1}\dots {j}_{d}\dots {i}_{D}}^{d}\left(\tau ,{\tau }^{\prime }\right)\right\vert }_{\text{SP}}\mathrm{d}\tau -1.\end{equation} \tag{ C10 }$$

The functional derivative Ξ^d i (τ, τ') can be obtained by varying the equation of motion of ${\xi }_{\boldsymbol{i}}^{+}$, as in the one-dimensional case:

$$\begin{equation}{{\Xi}}_{\boldsymbol{i}}^{d}\left(\tau ,{\tau }^{\prime }\right)={J}_{d}{\xi }_{\boldsymbol{i}}^{+}\left({\tau }^{\prime }\right)\theta \left(\tau -{\tau }^{\prime }\right){\mathrm{e}}^{{\int }_{{\tau }^{\prime }}^{\tau }\mathrm{d}s\left[-{\Gamma}{\xi }_{\boldsymbol{i}}^{+}\left(s\right)+\sum _{d=1}^{D}{J}_{d}{\varphi }_{\boldsymbol{i}}^{d}\left(s\right)\right]}.\end{equation} \tag{ C11 }$$

For a translationally invariant system one has ${\xi }_{\boldsymbol{i}}{\vert }_{\text{SP}}={\xi }_{\text{SP}}^{+}$, ${{\Xi}}_{\boldsymbol{i},\boldsymbol{j}}^{d}{\vert }_{\text{SP}}={{\Xi}}_{\text{SP}}^{d}$, ${\varphi }_{\boldsymbol{i}}^{d}{\vert }_{\text{SP}}={\varphi }_{\text{SP}}^{d}$, such that the SP equation simplifies to

$$\begin{equation}{\varphi }_{\text{SP}}^{d}\left({\tau }^{\prime }\right)=\frac{{\Gamma}}{{J}_{d}}{\int }_{0}^{{\tau }_{f}}{\left.{{\Xi}}_{\text{SP}}^{d}\left(\tau ,{\tau }^{\prime }\right)\right\vert }_{\text{SP}}\mathrm{d}\tau -1.\end{equation} \tag{ C12 }$$

For a fully isotropic system with J₁ =⋯= J_D = J, one additionally has ${\varphi }_{\text{SP}}^{d}={\varphi }_{\text{SP}}$ and the SP equations further simplify to

$$\begin{equation}{\varphi }_{\text{SP}}\left({\tau }^{\prime }\right)=\frac{{\Gamma}}{J}{\int }_{0}^{{\tau }_{f}}{\left.{{\Xi}}_{\text{SP}}\left(\tau ,{\tau }^{\prime }\right)\right\vert }_{\text{SP}}\mathrm{d}\tau -1,\end{equation} \tag{ C13 }$$

$$\begin{equation}{{\Xi}}_{\text{SP}}\left(\tau ,{\tau }^{\prime }\right)=\theta \left(\tau -{\tau }^{\prime }\right){\xi }_{\text{SP}}^{+}\left({\tau }^{\prime }\right){\mathrm{e}}^{{\int }_{{\tau }^{\prime }}^{\tau }\mathrm{d}s\left[-{\Gamma}{\xi }_{\text{SP}}^{+}\left(s\right)+DJ{\varphi }_{\text{SP}}\left(s\right)\right]},\end{equation} \tag{ C14 }$$

$$\begin{equation}{\dot {\xi }}_{\text{SP}}^{+}=\frac{{\Gamma}}{2}\left(1-{\xi }_{\text{SP}}^{+2}\right)+DJ{\xi }_{\text{SP}}^{+}{\varphi }_{\text{SP}}.\end{equation} \tag{ C15 }$$

Plateau equations can be derived from equations (C13) and (C15) in the large τ_f limit:

$$\begin{equation}{\varphi }_{\mathrm{P}}=\frac{DJ{\varphi }_{\mathrm{P}}}{{\Gamma}{\xi }_{\mathrm{P}}^{+}-DJ{\varphi }_{\mathrm{P}}},\end{equation} \tag{ C16a }$$

$$\begin{equation}{\varphi }_{\mathrm{P}}=-\frac{{\Gamma}}{2DJ}\frac{1-{\xi }_{\mathrm{P}}^{+2}}{{\xi }_{\mathrm{P}}^{+}}.\end{equation} \tag{ C16b }$$

These equations are solved by (34).

Here we provide a toy example illustrating the expansion of an integral which has two different SPs. We consider the integral

$$\begin{equation}I\left(a\right)={\int }_{-\infty }^{\infty }f\left(x\right)\mathrm{d}x,\end{equation} \tag{ C17 }$$

$$\begin{equation}f\left(x\right)={\mathcal{C}}_{a}\enspace {\mathrm{e}}^{-S\left(x\right)},\end{equation} \tag{ C18 }$$

$$\begin{equation}S\left(x\right)=\frac{{x}^{4}}{4{a}^{2}}-\frac{{x}^{2}}{2},\end{equation} \tag{ C19 }$$

where ${\mathcal{C}}_{a}$ is a normalization constant defined by I(a) = 1. Extremization of S(x) with respect to x yields two minima, x_SP = ±a. One can then expand the action around each SP as

$$\begin{equation}S={S}_{\text{SP}}+\frac{1}{2!}{S}_{\text{SP}}^{\left(2\right)}{\left(x-{x}_{\text{SP}}\right)}^{2}+{S}^{h},\end{equation} \tag{ C20 }$$

where S⁽²⁾ is the second variation evaluated at the SP and S^h includes all contributions of higher order. We then approximate equation (C17) as

$$\begin{equation}I\left(a\right)\approx {\mathcal{C}}_{a}\sum _{\text{s.p.}}{\mathrm{e}}^{-{S}_{\text{SP}}}{\int }_{-\infty }^{\infty }{\mathrm{e}}^{-\frac{1}{2}{S}_{\text{SP}}^{\left(2\right)}{\left(x-{x}_{\text{SP}}\right)}^{2}}\left[1+\sum _{m=3}^{n}{\alpha }_{m}{\left(x-{x}_{\text{SP}}\right)}^{m}\right]\mathrm{d}x,\end{equation} \tag{ C21 }$$

where the coefficients α_m are obtained by Taylor expanding ${\mathrm{e}}^{{S}^{h}}$ and the leftmost sum runs over the two SPs. For n ⩽ 2, none of the α_m is included and thus equation (C21) reduces to the evaluation of Gaussian fluctuations around the SP. To show how well the approximation (C21) captures the true value of I(a), in figure 11 we compare the exact and approximate integrands for different values of a, n. We find that the approximation (C21) for fixed n gets more accurate as a increases, such that the SPs are better spaced out. This is an example of the 'small overlap' condition discussed in the main text: one can separately expand about the two SPs and add up the individual contributions of the expansions, provided that regions (in this case, along the x axis) which contribute significantly to one integral give negligible contribution to the other.

Zoom In Zoom Out Reset image size

In section 3.2, we provided a general discussion of how QPTs emerge from the field theoretical application of the disentanglement approach. Namely, due to equation (43), a quantum critical point corresponds to the value for which there is an abrupt change in the dominant contribution to the grand state energy. This can be visualized for the quantum Ising chain, for which an exact analytical solution is available [65]. In figure 12, we show the small-Γ and large-Γ perturbative expansions of the ground state energy (44), given by (45). For any finite order in perturbation theory, the series cross at three points. As the order of both expansions is increased, the three crossing points converge toward a single point, the critical point Γ_c = J/2; see figure 12(a). Figure 12(b) further highlights that the small-Γ expansion is divergent for Γ > Γ_c; therefore, the relative terms will not contribute in this regime. Similarly, the large-Γ series does not contribute when Γ < Γ_c.

Zoom In Zoom Out Reset image size