Schwinger-Keldysh on the lattice: a faster algorithm and its application to field theory

TL;DR

This paper introduces a faster Monte Carlo algorithm for real-time field theory simulations that reduces computational costs and improves proposal efficiency, demonstrated on a 1+1 dimensional scalar ^4 theory.

Contribution

A novel algorithm that avoids explicit Jacobian calculations, significantly reducing computational complexity in real-time lattice field theory simulations.

Findings

Reduced update cost from O(N^3) to less than O(N^2)

Improved Monte Carlo proposal efficiency

Successfully applied to scalar ^4 theory with various couplings

Abstract

A new algorithm is developed allowing the Monte Carlo study of a 1 + 1 dimensional theory in real time. The main algorithmic development is to avoid the explicit calculation of the Jacobian matrix and its determinant in the update process. This improvement has a wide applicability and reduces the cost of the update in thimble-inspired calculations from O(N^3) to less than O(N^2). As an additional feature, the algorithm leads to improved Monte Carlo proposals. We exemplify the use of the algorithm to the real time dynamics of a scalar {\phi}^4 theory with weak and strong couplings.

Tables

Table 1. Table 1: Simulation parameters for the ensembles used in this study. The last two sets of columns indicate the number of updates required to generate the ensembles used in this study, the step-size δ 𝛿 \delta , and the acceptance rate. Grady indicates the algorithm proposed in this paper, and J 0 subscript 𝐽 0 J_{0} indicates the variant of the old algorithm described at the end of Section  IV .

$\lambda$		$T_{\text{flow}}$	$\mathrm{Re}{⟨\Phi ⟩}_P$	$J_0$				Grady
				$\mathrm{updates}/{10}^6$	$\delta$	acceptance	stat power	$\text{updates}/{10}^6$	$\delta$	acceptance
0.1		1.8	0.278(2)	100	0.11	53%	0.99	3.2	0.1	57%
0.5		1.6	0.193(2)	50	0.115	50%	0.90	2.4	0.1	57%
1.0		1.6	0.189(2)	50	0.1	53%	0.68	2.4	0.1	55%

Equations

$$⟨\hat{\mathcal{O}}(t)\hat{\mathcal{O}}(t^{\prime })⟩=\mathrm{Tr}\left(\hat{\mathcal{O}}(t)\hat{\mathcal{O}}(t^{\prime })\hat{\rho }\right).$$

$$⟨{\hat{\mathcal{O}}}_1(x_1){\hat{\mathcal{O}}}_2(x_2)\cdots ⟩=\mathrm{Tr}[\hat{\rho }(0){\hat{\mathcal{O}}}_1(x_1){\hat{\mathcal{O}}}_2(x_2)\cdots ],$$

$$Z[J_{+},J_{-},J_{\beta }]=\mathrm{Tr}[U(T-i\beta ,T;J_{\beta })U(T,T^{\prime };J_{-})U(T^{\prime },T;J_{+})],$$

$$Z[J_{+},J_{-},J_{\beta }]=\int D{\varphi }_{+}D{\varphi }_{-}D{\varphi }_{\beta }{e}^{i{S}_{SK}[{\varphi }_{+},{\varphi }_{-},{\varphi }_{\beta },J_{+},J_{-},J_{\beta }]}$$

$$S_{SK}[{\varphi }_{+},{\varphi }_{-},{\varphi }_{\beta },J_{+},J_{-},J_{\beta }]$$

$$⟨T\varphi (t_1,{\mathbf{x}}_1)\varphi (t_2,{\mathbf{x}}_2)⟩=\frac{{\delta }^{2}Z[J_{+},J_{-},J_{\beta }]}{\delta J_{+}(t_1,{\mathbf{x}}_1)\delta J_{+}(t_2,{\mathbf{x}}_2)}.$$

$$\begin{array}{rl}S[\varphi ]\equiv -i{S}_{SK,lattice}=\sum _{t,n}{a}_{t}a& [\frac{({\varphi }_{t+1,n}-{\varphi }_{t,n}{)}^2}{2a_t^2}+\frac{1}{2}\left(\frac{({\varphi }_{t+1,n+1}-{\varphi }_{t+1,n}{)}^2}{2a^2}+\frac{({\varphi }_{t,n+1}-{\varphi }_{t,n}{)}^2}{2a^2}\right)\\ & +\frac{1}{2}{m}^2\frac{{\varphi }_{t,n}^2+{\varphi }_{t+1,n}^2}{2}+\frac{\lambda }{4!}\frac{{\varphi }_{t+1,n}^4+{\varphi }_{t,n}^4}{2}],\end{array}$$

$$\begin{array}{rl}{a}_t=& ia,\text{for}0\le t<N_t,\\ a_t=& a,\text{for}{N}_t\le t<N_t+N_{\beta }/2,\\ a_t=& -ia,\text{for}{N}_t+N_{\beta }/2\le t<2N_t+N_{\beta }/2,\\ a_t=& a,\text{for}2N_t+N_{\beta }/2\le t<2N_t+N_{\beta },\end{array}$$

$$\displaystyle\langle\phi_{t_{1},n_{1}}\phi_{t_{2},n_{2}}\rangle={\int\big{(}\prod_{t,n}d\phi_{t,n}\big{)}e^{-S[\phi]}\phi_{t_{1},n_{1}}\phi_{t_{2},n_{2}}\over\int\big{(}\prod_{t,n}d\phi_{t,n}\big{)}e^{-S[\phi]}}\,.$$

$$\begin{array}{rl}⟨{\varphi }_i{\varphi }_j⟩& \approx ⟨{\varphi }_i{\varphi }_j{⟩}_0+\lambda \left(-\frac{1}{4!}\right)\sum _k{\tilde{a}}_k⟨{\varphi }_i{\varphi }_j{\varphi }_k^4{⟩}_0^c+{\lambda }^2\frac{1}{2}{\left(-\frac{1}{4!}\right)}^2\sum _{k,l}{\tilde{a}}_k{\tilde{a}}_l⟨{\varphi }_i{\varphi }_j{\varphi }_k^4{\varphi }_l^4{⟩}_0^c\\ & =(H_0{)}_{ij}^{-1}-\frac{\lambda }{2}\sum _k{\tilde{a}}_k(H_0{)}_{ik}^{-1}(H_0{)}_{kk}^{-1}(H_0{)}_{kj}^{-1}+\frac{{\lambda }^2}{4}\sum _{k,l}{\tilde{a}}_k{\tilde{a}}_l[(H_0{)}_{il}^{-1}((H_0{)}_{lk}^{-1}{)}^2(H_0{)}_{kk}^{-1}(H_0{)}_{kj}^{-1}\\ & +(H_0{)}_{il}^{-1}(H_0{)}_{lk}^{-1}(H_0{)}_{ll}^{-1}(H_0{)}_{kk}^{-1}(H_0{)}_{kj}^{-1}+\frac{2}{3}(H_0{)}_{il}^{-1}((H_0{)}_{lk}^{-1}{)}^3(H_0{)}_{kj}^{-1}],\end{array}$$

$$\frac{d{z}^i(\tau )}{d\tau }=\overline{\frac{\partial S(z(\tau ))}{\partial z^i}}\text{with}z(0)=x\in {\text{mathds}R}^N.$$

$$\frac{d{\omega }^i(\tau )}{d\tau }=\overline{H_{ij}(z(\tau )){\omega }^j(\tau )}\text{with}\omega (0)=v\text{and}{H}_{ij}\equiv \frac{{\partial }^{2}S}{\partial z^i\partial z^j}.$$

$$J_{ij}\equiv \mathrm{T}f(e_j{)}^i,$$

$$\frac{d}{d\tau }{⟨w,v⟩}_{\text{mathds}C}=\frac{d}{d\tau }({\overline{w}}^i{v}^i)={\overline{\frac{dw}{d\tau }}}^i{v}^i+\overline{w^i}{\frac{dv}{d\tau }}^i=H_{ij}{w}^j{v}^i+\overline{w^i{H}_{ij}{v}^j}=2\text{Re}[w^i{H}_{ij}{v}^j].$$

$$\mathrm{Im}{⟨\mathrm{T}f(e_i),\mathrm{T}f(e_j)⟩}_{\text{mathds}C}=\mathrm{Im}{⟨e_i,e_j⟩}_{\text{mathds}C}=0\text{which can be restated as}\mathrm{Im}\overline{J_{ki}}{J}_{kj}=\mathrm{Im}(J^{†}J{)}_{ij}=0.$$

$${⟨\mathrm{T}f(e_j),i\mathrm{T}f(e_k)⟩}_{\text{mathds}R}=\mathrm{Re}{⟨\mathrm{T}f(e_j),i\mathrm{T}f(e_k)⟩}_{\text{mathds}C}=-\mathrm{Im}{⟨\mathrm{T}f(e_j),\mathrm{T}f(e_k)⟩}_{\text{mathds}C}=0.$$

$$v^{\prime }=v_{\parallel }^{\prime }+v_{\perp }^{\prime }\text{with}{v}_{\parallel }^{\prime }=v_{\parallel }^j\mathrm{T}f(e_j)=\mathrm{T}f(v_{\parallel }^j{e}_j)\text{and}{v}_{\perp }^{\prime }=v_{\perp }^{j}i\mathrm{T}f(e_j)=i\mathrm{T}f(v_{\perp }^j{e}_j),$$

$$Z={\int }_{\mathcal{M}}\mathrm{d}{x}^{\prime }{e}^{-S(x^{\prime })}={\int }_{{\text{mathds}R}^N}\mathrm{d}x\det J(x)e^{-S(f(x))}={\int }_{{\text{mathds}R}^N}\mathrm{d}x|\det J(x)|e^{-\mathrm{Re}S(f(x))}\Phi (x),$$

$$⟨\mathcal{O}(x)⟩=\frac{{⟨\mathcal{O}(x)\Phi (x)⟩}_P}{{⟨\Phi (x)⟩}_P}\text{with}P(x)=\left|\det J(x)\right|e^{-\mathrm{Re}S(f(x))}.$$

$$\mathrm{Pr}(x_n\to x_{n+1})=\sqrt{\frac{\det (J^{†}J{)}_n}{{\pi }^N{\delta }^2}}{e}^{-{\Delta }^T(J^{†}J{)}_n\Delta /{\delta }^2},$$

$$P(\xi )=\frac{\det (J^{†}J{)}_{n+1}}{{\pi }^N}{e}^{-{\xi }^{†}(J^{†}J{)}_{n+1}\xi },$$

$${\mathrm{P}}_{\text{acc}}=\min \{1,e^{-\Delta S+{\xi }^{†}[(J^{†}J{)}_{n+1}-(J^{†}J{)}_n]\xi +{\Delta }^T[(J^{†}J{)}_n-(J^{†}J{)}_{n+1}]\Delta /{\delta }^2}\}.$$

$$T(x_n\to x_{n+1})=\mathrm{Pr}(x_n\to x_{n+1})\int \mathrm{d}\xi \mathrm{d}{\xi }^{†}P(\xi ){\mathrm{P}}_{\text{acc}}.$$

$$\begin{array}{rl}P(x_n)T& (x_n\to x_{n+1})=e^{-\mathrm{Re}S(x_n^{\prime })}\left|\det J_n\right|\sqrt{\frac{\det (J^{†}J{)}_n}{{\pi }^N{\delta }^2}}{e}^{-{\Delta }^T(J^{†}J{)}_n\Delta /{\delta }^2}\\ & \times \int \mathrm{d}\xi \mathrm{d}{\xi }^{†}\frac{\det (J^{†}J{)}_{n+1}}{{\pi }^N}{e}^{-{\xi }^{†}(J^{†}J{)}_{n+1}\xi }\min \{1,e^{-\Delta S+{\xi }^{†}[(J^{†}J{)}_{n+1}-(J^{†}J{)}_n]\xi +{\Delta }^T[(J^{†}J{)}_n-(J^{†}J{)}_{n+1}]\Delta /{\delta }^2}\}\\ =& \frac{\det (J^{†}J{)}_n\det (J^{†}J{)}_{n+1}}{\delta {\pi }^{3N/2}}\\ & \times \int \mathrm{d}\xi \mathrm{d}{\xi }^{†}\min \{e^{-\mathrm{Re}S(x_n^{\prime })-{\Delta }^T(J^{†}J{)}_n\Delta /{\delta }^2-{\xi }^{†}(J^{†}J{)}_{n+1}\xi },e^{-\mathrm{Re}S(x_{n+1}^{\prime })-{\Delta }^T(J^{†}J{)}_{n+1}\Delta /{\delta }^2-{\xi }^{†}(J^{†}J{)}_n\xi }\}.\end{array}$$

$$⟨\mathcal{O}⟩=\frac{{⟨\mathcal{O}(x)\tilde{\Phi }(x)⟩}_{P_0}}{{⟨\tilde{\Phi }(x)⟩}_{P_0}}\text{with}\tilde{\Phi }(x)=\frac{\det J(x)}{\left|\det J(0)\right|}{e}^{-i\mathrm{Im}S(f(x))},$$

$$C(t,p)={⟨\varphi (t,p)\varphi (0,p{)}^{†}⟩}_{\beta }\text{with}\varphi (t,p)\equiv \frac{1}{N_x}\sum _{n=0}^{N_x-1}{e}^{ipn}{\varphi }_{t,n}.$$

$$C_b(t,p)={⟨\varphi (0,p{)}^{†}\varphi (t,p)⟩}_{\beta },$$

$$C(t,p)=\frac{1}{2(N_t+1-t)}\sum _{t^{\prime }=0}^{N_t-t}\left[⟨\varphi (t+t^{\prime },p)\varphi (t^{\prime },-p)⟩+\overline{⟨\varphi (t^{\prime }+T_b,p)\varphi (t+t^{\prime }+T_b,-p)⟩}\right]\text{ for }t=0,\dots ,N_t,$$

$$\text{stat power}=\frac{1}{N_{\text{cfg}}}\frac{{\left({\sum }_i|\det J({\varphi }_i)|\right)}^2}{{\sum }_i|\det J({\varphi }_i){|}^2},$$

Full text

I Introduction

Some of the most interesting – and challenging – problems in many-body physics are dynamical questions. They describe near equilibrium states through transport coefficients like diffusion constants, viscosities, conductivities, as well as phenomena away from equilibrium. The central objects of interest in these class of problems are the time dependent correlation functions of the form

[TABLE]

where the operators have the usual time evolution $\mathcal{\hat{O}}(t)=e^{i\hat{H}t}\mathcal{\hat{O}}(0)e^{-i\hat{H}t}$ and $\hat{\rho}$ is the density matrix which reduces to the Boltzmann form, $e^{-\beta\hat{H}}/{\rm Tr}(e^{-\beta\hat{H}})$, in equilibrium. Problems in this category show up in almost every field of physics: cosmology, heavy ions collisions, and condensed matter physics to name a few.

Unfortunately the tools for tackling such dynamical problems from first principles are very limited. Even in weakly coupled systems, the study of long time dynamics (or low momentum properties) require complicated resummations of the perturbative expansion Braaten and Pisarski (1990); Jeon and Yaffe (1996). Monte-Carlo techniques, the method of choice for non-perturbative problems, have a fundamental difficulty in dealing with real time (as opposed to imaginary time) dynamics, due to a particularly severe version of the “sign problem”. Many Monte Carlo-based approaches, including the one used in relativistic theories and that we use in this paper, are based on a path integral representation of the observable of interest. Such a path integral representation exists for real time observables in or out of equilibrium, and is based on the the Schwinger-Keldysh formalism Schwinger (1961); Keldysh (1964). The problem in the Monte-Carlo evaluation of this path integral is that the integrand is a pure phase, as opposed to a fast decaying real function, and the importance sampling of the integrand, based on the interpretation of the integrand as a probability density, is not possible. As we will comment below, the sign problem for path integrals in the Schwinger-Keldysh formalism is, in a certain sense, the worst possible.

We are aware of two approaches that address this problem via Monte-Carlo techniques. The first one is to concentrate on near equilibrium and attempt to compute transport coefficients. They can be computed through the Kubo formula from the knowledge of certain equilibrium real time correlators (1). In principle, the correlators in imaginary time, $\langle\mathcal{\hat{O}}_{E}(\tau)\mathcal{\hat{O}}_{E}(\tau^{\prime})\rangle$ with $\mathcal{\hat{O}}_{E}(\tau)=e^{\hat{H}\tau}\mathcal{\hat{O}}(0)e^{-\hat{H}\tau}$, contain the same information as the real time ones and can be computed with standard Monte Carlo techniques, frequently without a sign problem Karsch and Wyld (1987); Huang (1993); Meyer (2008, 2007); Ding et al. (2011). In practice, however, exponentially good precision on imaginary time is required to reconstruct it on real time. The second approach is to use Langevin methods (“stochastic quantization”Parisi and s. Wu (1981)). The drawback of the (complex) Langevin approach is that it does not always converge, or sometimes converges to an incorrect result 111For a recently proposed criterion for convergence see Nagata et al. (2016). In fact, the few attempts of applying the complex Langevin method to real time dynamics seem to suggest that it converges to the wrong result if the time separation $t-t^{\prime}$ is more than the inverse temperature $\beta$ Berges and Stamatescu (2005); Berges et al. (2007); Mizutani et al. (2008).

In the last few years a new approach to compute path integrals with sign problem has been developed Cristoforetti et al. (2012). Although different versions vary in detail, they are all based in the deformation of the path integral from real values of the fields to a suitably chosen middle dimensional (i.e. with the same dimensions as the real field space) submanifold of the complexified field space. The equality of the integral over this new manifold to the integral over the original real space is guaranteed by a multidimensional version of Cauchy’s theorem and by choosing the asymptotic properties of the manifold properly (the analogue of avoiding the “arcs at infinity” familiar from complex analysis of functions of one complex variable). One choice of such a manifold is to deform the contour of integration from the real space to an appropriate combination of “thimbles”, the multidimensional analogues of the “steepest descent” or “constant phase” path from the theory of one complex variable. The sign problem is solved because, along thimbles, the phase of the integrand is constant. A difficulty with this choice of manifold is that it is, in general, nearly impossible to determine the particular combination of thimbles that is equivalent to the original region of integration. There are also the issues of how to sample disconnected thimbles and find them as their location is not known analytically.

Another choice of integration manifold was proposed in Alexandru et al. (2016a) and pursued by our group recently Alexandru et al. (2016b, 2017a, 2017b). In this proposal, the manifold of integration is obtained by taking the real fields as a starting value and evolving them according to the (anti)holomorphic gradient flow (the complex conjugate of the gradient of the action). This flow evolves a given real field configuration along a particular trajectory, determined by the gradient of the action, in the complexified field space. The end point of this evolution is determined by the “flow time”, which is viewed as a free parameter. Therefore flowing the original real field space by some flow time creates an alternative, complex path integration manifold associated with the value of that flow time. In the limit of large flow times this manifold coincides with the precise combination of thimbles which is equivalent to the original integration domain. For finite flow times it provides a manifold i) that is equivalent to the original domain of integration, ii) on which the phase variation of the integrand is milder than on the real space and iii) that is connected, making the stochastic sampling easier to accomplish.

This method was applied to the equilibrium real time dynamics of an anharmonic oscillator in Alexandru et al. (2016b). The correct result (which is known in this case through direct diagonalization of the hamiltonian) was obtained, even for time differences $t-t^{\prime}$ of the order of $\approx 4\beta$. However, the particular implementation of this method had two major shortcomings which prevented a similar computation in a field theory. Firstly, an unreasonably large number of Monte-Carlo steps were required for thermalization and decorrelation. The reason for this was traced back to the fact that in this method, it is natural to make isotropic proposals in real space that are then “flowed” to the manifold of integration where they are highly anisotropic. In Alexandru et al. (2016b) some attempts were made to correct for this anisotropy by using a combination of gaussian approximation and trial-and-error adjustment for the proposals to make them more isotropic when flowed to the manifold of integration, but with limited success. It also required the computation and storage of a set of $N$, $N$-dimensional vectors (where $N$ is the number of degrees of freedom on the lattice), which was possible in the anharmonic oscillator problem but is prohibitive in a field theory with a large lattice. Secondly, the need to compute the jacobian associated with the parametrization of the manifold by its real coordinates makes every step of the Monte-Carlo chain computationally expensive. Previously, this problem had been dealt with by using an estimator of the jacobian Alexandru et al. (2016c, d, 2017a) and reweighting the difference when making measurements. The estimator we developed is likely to be useful if the coupling is small and/or the manifold of integrations is nearly parallel to the real plane. This is not the case for the real time calculations so a new method that bypasses the need to compute the jacobian at every step of the Monte Carlo chain is necessary.

The purpose of this paper is to present an algorithm without these two difficulties, and one that does not require the storage of the $N$, $N$-dimensional vectors. The main idea is the Grady algorithm Grady (1985); Creutz (1992) which is used in lattice QCD in order to avoid the computation of fermion determinants. The effect of the jacobian is embedded in a bias of the proposals that are isotropic in the flowed manifold. We implement the Grady algorithm into the holomorphic gradient flow method to perform a real time calculation on a $\phi^{4}$ theory in $1+1$ dimensions. In this model and for weak or strong couplings, we are able to perform an even computationally cheaper calculation by approximating the holomorphic flow by its gaussian approximation in the calculation of the proposal (and reweighting the difference between them when making measurements).

In Section 2 we briefly review properties of the holomorphic flow and the new algorithm is described. In Section 3 the model and its representation in the Schwinger-Keldysh formalism is described. Results are described in Section 4 and a conclusion is presented in Section 5.

II The Schwinger-Keldysh path integral and the model

As discussed in the introduction, we are interested in expectation values of the form

[TABLE]

where $\mathcal{\hat{O}}_{i}(x)=e^{i\hat{H}t-i{\bf p}\cdot{\bf x}}\mathcal{\hat{O}}_{i}(t=0,{\bf x}={\bf 0})e^{-i\hat{H}t+i{\bf p}\cdot{\bf x}}$ are operators in the Heisenberg representation and $\hat{\rho}(0)=e^{-\beta\hat{H}(0)}/{{\rm Tr}(e^{-\beta\hat{H}(0)})}$ is the density matrix representing the initial state of the system. When the time evolution is determined by a time independent hamiltonian $\hat{H}(t)=\hat{H}(0)=\hat{H}$, the system is in thermodynamic equilibrium which we will assume. The non-equilibrium case can also be studied within the formalism after a slight modification. We will briefly comment on this point later. Expectation values of this form can be obtained from the generating functional

[TABLE]

where $J_{\pm},J_{\beta}$ are external classical currents coupled to the fields in the theory and $U(T,T^{\prime};J)$ is the time evolution operator under the influence of the external current $J$ between times $T$ and $T^{\prime}$. In order to compute correlators as in (2) we need $T<t_{1},t_{2},\cdots<T^{\prime}$. The generating function has the path integral representation Schwinger (1961); Keldysh (1964)

[TABLE]

where the action, $S_{SK}$, is defined as an integral of the Lagrangian along a time contour, ${\cal C}$, which lives in the complex plane:

[TABLE]

with the boundary conditions $\phi_{+}(T^{\prime})=\phi_{\beta}(T^{\prime}),\phi_{\beta}(T^{\prime}-i\beta/2)=\phi_{-}(T^{\prime}-i\beta/2),\phi_{-}(T-i\beta/2)=\phi_{\beta}(T-i\beta/2),\phi_{\beta}(T-i\beta)=\phi_{+}(T)$ ($\phi_{\beta}(T-i\beta)=-\phi_{+}(T)$ for fermionic fields).

The contour ${\cal C}$ is depicted in Fig. 1 (left). The upper and lower parts of the contour, parallel to the real axis, are associated with the forward and backward time evolution (i.e. the second and third terms in the trace in Eq. 3). The parts that are along the imaginary axis are associated with the insertion of the density matrix (i.e. the first term in the trace in Eq. 3). Note that we chose to split the density matrix into two parts that are inserted at times $T$ and $T^{\prime}$. This choice assumes the existence of equilibrium where the Hamiltonian is time independent and the density matrix commutes with time evolution operator. In order to study an out-of-equlibrium system this contour has to be modified such that the density matrix is inserted at time $T^{\prime}$ as a whole. Even though our construction can be generalized, we will not discuss this case in this paper.

Given the path integral representation, the various correlators can be computed by differentiating $Z[J_{+},J_{-},J_{\beta}]$ with respect to external sources. The time ordering is such that the operators that are inserted in the lower branch of ${\cal C}$ always have a larger time compared to those that are inserted at the upper branch. A two point correlator with both operators are in the upper (lower) branch is time ordered (anti time ordered). For instance

[TABLE]

In the paper, we present a Monte Carlo method to compute time dependent correlation functions as in (7) for the $1+1$ dimensional $\phi^{4}$ theory with the potential $V[\phi]={1\over 2}m^{2}\phi^{2}+{\lambda\over 4!}\phi^{4}$. The lattice action corresponding to the Schwinger-Keldysh path integral is given by

[TABLE]

where $t$ and $n$ indexes the lattice along the time and spatial directions, $a$ is the spatial lattice spacing and $a_{t}$ is the time lattice spacing:

[TABLE]

and $N_{t},N_{\beta}$ are the number of lattice points on the real and imaginary axis, respectively. This is shown in Fig. 1 (right). We embedded a factor of $-i$ in the definition of the action in Eq. 8 so that the measure in the path integral is $e^{-S}$. The change in $a_{t}$ with $t$ determines the contour in the complex time plane defining the Schwinger-Keldysh action. The fields $\phi_{+},\phi_{-}$ and $\phi_{\beta}$ correspond, in our discretized action, to $\phi_{t,n}$ for $t$ in the ranges $0\leq t<N_{t}$, $N_{t}+N_{\beta}/2\leq t<2N_{t}+N_{\beta}/2$, and $N_{t}\leq t<N_{t}+N_{\beta}/2$ and $2N_{t}+N_{\beta}/2\leq t<2N_{t}+N_{\beta}/2$ respectively. The correlators that we are interested in computing are now given by the discretized path integral

[TABLE]

Along the two branches of the real axis the measure, $e^{-S[\phi]}$, of the path integral is a pure phase, and lacks exponential damping as the value of $\phi_{\pm}(t,\mathbf{x})$ is varied, in contrast to the Euclidean branch. Thus, any attempt at reweighting the phase of the integrand is doomed to fail as the average phase vanishes. In this sense, the sign problem arising in real time is the worst possible. In the following two sections, we will present a new Monte-Carlo algorithm to alleviate this sign problem by using the holomorphic gradient flow.

In order to validate our results at small coupling and to ascertain that our results at larger coupling indeed cannot be reproduced by perturbation theory, we compare against perturbative results. The propagator in the Schwinger-Keldysh contour is given up to order ${\cal O}(\lambda^{2})$ by

[TABLE]

where $i,j,k,l$ are combined time and space indices, $\tilde{a}_{k}=a(a_{t(k)-1}+a_{t(k)})/2$, $H_{0}$ is the Hessian of the Schwinger-Keldysh action at $\phi=0$ and $\langle\cdot\rangle_{0}^{c}$ denotes the connected part of the correlator.

III Holomorphic gradient flow

In order to solve the sign problem, that is to reduce the phase fluctuations for the path integral to a level where we can easily reweight it, we will deform the integration manifold of our path integral. The first step of the process is to promote all real variables (values of the field at each space-time point) to complex ones. Using a generalized version of Cauchy’s theorem, we can show that we can deform the integration manifold in the complex space without changing the value of the integral, as long as we do not cross any singularities of the integrand and we preserve the asymptotic behavior of the field. There is a rather large latitude in choosing the deformation but in this work we will use a deformation induced by the holomorphic gradient flow. In this section we will review briefly the relevant details.

The holomorphic gradient flow is defined, for a system with $N$ real degrees of freedom, through the set of differential equations

[TABLE]

Integrating the flow equation above for a fixed amount of “time” $T_{\text{flow}}$ defines a map $x\to f(x)$ where $f(x)\equiv z(T_{\text{flow}})$. The image under this map of the original integration domain ${\mathds{R}}^{N}$ is our new integration manifold ${\cal M}=f({\mathds{R}}^{N})$. Note that this manifold depends on $T_{\text{flow}}$ and as the flow time is increased the sign fluctuations become milder Alexandru et al. (2016e). When the action is real the flow keeps the points in the real subspace, but for complex actions the image manifold ${\cal M}$ will be different from ${\mathds{R}}^{N}$. Since we use the points in ${\mathds{R}}^{N}$ to keep track of the points in the integration manifold ${\cal M}$, we will refer to ${\mathds{R}}^{N}$ as the parametrization manifold.

The map $f$ naturally induces a linear map, denoted by ${\rm T}\!f$, between the tangent space at point $x$ and the tangent space at point $x^{\prime}=f(x)$. A vector $v$ tangent at $x$ is mapped to ${\rm T}\!f(v)=v^{\prime}$ tangent at $x^{\prime}$ with $v^{\prime}=\omega(T_{\text{flow}})$ where $\omega$ satisfies the differential equation

[TABLE]

This equation is derived by considering the flow of a point infinitesimally displaced in direction $v$ away from $x$. Note that in the equation above the Hessian $H$ is to be evaluated along the trajectory $z(\tau)$ that takes $x$ to $x^{\prime}$, so this equation needs to be solved together with Eq. 12 in order to “transport” a vector.

Since the vector flow is a linear map, it can be represented by a matrix. Taking in ${\mathds{R}}^{N}$ the canonical basis $\{e_{j}\}_{j=1,\ldots,N}$ with $e_{j}^{i}=\delta_{ij}$ the vectors ${\rm T}\!f(e_{j})$ will form a basis for the tangent space of $\mathcal{M}$ at $x^{\prime}$. The matrix $J$ which has these vectors as columns, that is

[TABLE]

can be used to map $v=v^{i}e_{i}$ into $v^{\prime}=v^{i}{\rm T}\!f(e_{i})$ so that the components of $v^{\prime}$ are $v^{\prime i}=J_{ij}v^{j}$. Note that while the components of $v$ have to be real, the components of $v^{\prime}$ are in general complex, since the tangent space at $x^{\prime}$ is usually not parallel to ${\mathds{R}}^{N}$. We stress that $v^{\prime i}$ are the components of $v^{\prime}$ in the canonical basis $e_{i}$ in ${\mathds{C}}^{N}$, not the components of $v^{\prime}$ in the basis ${\rm T}\!f(e_{i})$. If we decompose $v^{\prime}$ in the basis ${\rm T}\!f(e_{i})$ the components are real, as they should be, and they are in fact $v^{i}$. The map ${\rm T}\!f$ and the matrix $J$ depend on the starting point $x\in{\mathds{R}}^{n}$, but to simplify the notation we will make this explicit only when required.

An important property of the tangent map ${\rm T}\!f$ is that it preserves the imaginary part of the dot product of two vectors since the derivative of the product along the flow is real:

[TABLE]

We denoted the ${\mathds{C}}^{N}$ scalar product with $\left\langle w,v\right\rangle_{\mathds{C}}$ to distinguish it from the real scalar product $\left\langle w,v\right\rangle_{\mathds{R}}=\operatorname{Re}w^{i}\operatorname{Re}v^{i}+\operatorname{Im}w^{i}\operatorname{Im}v^{i}=\operatorname{Re}\left\langle w,v\right\rangle_{\mathds{C}}$. The invariance of the imaginary part implies that for any $i$ and $j$

[TABLE]

Thus the matrix $J^{\dagger}(x)J(x)$ is real. Another important property is that while $\{{\rm T}\!f(e_{j})\}_{j=1,\ldots,N}$ span the tangent space at $x^{\prime}$ (viewed as a real vector space), the vectors $\{i{\rm T}\!f(e_{j})\}_{j=1,\ldots,N}$ span the orthogonal vector space. This can be seen by considering the real scalar product between any two vectors in these sets:

[TABLE]

This means that any vector $v^{\prime}$ at $x^{\prime}$ can be decomposed as

[TABLE]

where $v_{\parallel,\perp}^{j}$ coefficients are real. Defining $v_{\parallel,\perp}\equiv v_{\parallel,\perp}^{j}e_{j}$, we have $v^{\prime}_{\parallel}=Jv_{\parallel}$ and $v^{\prime}_{\perp}=iJv_{\perp}$.

The partition function is evaluated by integrating over the manifold ${\cal M}$:

[TABLE]

where $\Phi(x)\equiv\exp[-i\operatorname{Im}S(f(x))+i\arg\det\,J(x)]$ is a pure phase. Note that the measure $dx^{\prime}$ is the Cauchy measure on the manifold ${\cal M}$ embedded in ${\mathds{C}}^{N}$ and that the change of variables from the integration manifold to the parametrization manifold ${\mathds{R}}^{N}$ gives rises to the Jacobian $\det\,J(x)$, which is the determinant of the matrix $J$ defined in Eq. 14. To evaluate observable averages with respect to $Z$, we will sample configurations according to the positive weight $P(x)$ and then evaluate observables by reweighting, that is, we compute averages from the formula

[TABLE]

The averages $\left\langle\cdot\right\rangle_{P}$ are taken with respect to the probability weight $P(x)$. We discuss how to sample this measure in the next section.

IV Algorithm

In this section we will present a method of sampling configurations according to the probability weight $P(x)$ defined in Eq. 20. The method described here is based on the Metropolis algorithm, where new configurations are proposed and an accept-reject step is used to ensure detailed balance. The algorithm we used previously  Alexandru et al. (2016e) suffered from two main flaws. The first was the substantial cost of computing the jacobian $J(x)$ at every step of the Markov chain. The second was that the proposals, straightforwardly, were chosen to be isotropic in the real variables $x$. As those variables parametrize the actual manifold of integration $\mathcal{M}$ through the very non-linear map $x^{\prime}=f(x)$, the resulting proposals were very anisotropic in $\mathcal{M}$. This distortion effect was partially, but not completely, compensated in  Alexandru et al. (2016e). The result was that very small proposal steps had to be chosen in order to achieve reasonable acceptance rates and made the algorithm perform poorly. The algorithm we discuss here improves on the one used in  Alexandru et al. (2016e) in both respects.

To better understand the advantages of this algorithm, let us first review the computationally intensive steps involved. The cost estimates below will concern typical bosonic systems, as the one considered in this paper. The most expensive part of the computation of the action of the flowed configuration $S(x^{\prime})=S(f(x))$, is the calculation of the flowed configurations $x^{\prime}$, which requires the integration of Eq. 12. For this, we use an adaptive Runge-Kutta integrator Cash and Karp (1990). In terms of scaling with the size of the system this is a ${\cal O}(N)$ calculation, if we assume that $T_{\text{flow}}$ remains fixed as we increase the system size. To compute $J(x)$ we need to integrate Eq. 13 for each vector in the basis. The Hessian is usually a sparse matrix, so each integration step can be implemented with complexity ${\cal O}(N)$, and the cost of computing ${\rm T}\!f(v)$ for some vector $v$ is ${\cal O}(N)$. Overall, the cost of computing $J(x)$ is then ${\cal O}(N^{2})$ and its determinant has a cost of order ${\cal O}(N^{3})$. For large systems the cost of the later steps quickly becomes dominant. For some systems estimators of $\left|\det\,J(x)\right|$ can be employed for sampling and the difference can be reweighted Alexandru et al. (2016c). However, these estimators do not work well for the system considered in this paper.

The second problem is to find proposals that move efficiently through the integration manifold ${\cal M}$. Proposals that are distributed isotropically in the parametrization space are in general mapped to a highly skewed distribution in the integration space, leading to an inefficient sampling of the manifold Alexandru et al. (2016e). The skewed distribution appears because the map ${\rm T}\!f$ scales very differently vectors that point in different directions, that is, the eigenvectors of $J(x)$ have eigenvalues of very different magnitudes. When the matrix $J(x)$ is relatively constant over the sampled region in the parametrization space, we can bias the proposals in the parametrization space such that their distribution in the integration manifold is relatively isotropic. To be specific, denote the current configuration by $x_{n}$ and the proposed one by $x_{n+1}$. When the proposal is selected with probability ${\rm Pr}(x_{n}\to x_{n+1})\propto\exp(-\Delta^{T}M\Delta)$, where $\Delta\equiv x_{n+1}-x_{n}$ and $M$ is a fixed real positive-definite matrix, the probability ${\rm Pr}$ is symmetric in $x_{n}$ and $x_{n+1}$, so the acceptance probability required for detailed balance is the usual ${{\rm P}_{\text{acc}}}=\min\{1,P(x_{n+1})/P(x_{n})\}$. When the parametrization manifold is tangent to a critical point $x_{\text{cr}}$, a good choice for matrix $M$ is $J(x_{\text{cr}})^{\dagger}J(x_{\text{cr}})$, since this bias can be effectively constructed using the “eigenvectors” and “eigenvalues” of the Hessian evaluated at the critical point Alexandru et al. (2016e, d). When flowing from the original integration manifold ${\mathds{R}}^{N}$, as we will do in the present paper, a possible choice for the matrix $M$ would be the quadratic approximation to the real part of the action $\operatorname{Re}S(f(x))\approx-x^{T}Mx$ Alexandru et al. (2016b). The problems with these methods are that i) the accept-reject step requires the calculation of $J(x)$ and its determinant, and ii) they are only effective when the matrix $J$ does not fluctuate too much over the sampled configurations. Ideally, we would like to make proposals that are isotropic around $x_{n}^{\prime}\in{\cal M}$, that is ${\rm Pr}(x_{n}\to x_{n+1})\propto\exp(-\Delta^{T}J^{\dagger}(x_{n})J(x_{n})\Delta)$. These proposals are isotropic in the tangent space at $x_{n}^{\prime}$ because $\eta^{\prime}_{\parallel}\equiv J(x_{n})\Delta$ is a random vector in this space distributed with probability $P(\eta_{\parallel}^{\prime})\propto\exp(-\eta_{\parallel}^{\prime\dagger}\eta_{\parallel}^{\prime})$. These proposals are not symmetric under the exchange $x_{n}\leftrightarrow x_{n+1}$, since the matrices $J_{n}\equiv J(x_{n})$ and $J_{n+1}\equiv J(x_{n+1})$ are different. To account for this asymmetry, the acceptance probability needs to be modified to ${{\rm P}_{\text{acc}}}=\min\{1,P(x_{n+1})\left|\det J_{n+1}\right|/P(x_{n})\left|\det J_{n}\right|\}$. This is still expensive, since the determinants appearing in $P(x)$ do not cancel in the acceptance ratio and we are still required to compute $\left|\det J\right|$. However, this suggests a way to arrange the proposals to cancel the determinants: if we could make proposals using the probability ${\rm Pr}(x_{n}\to x_{n+1})\propto\exp[-\Delta^{T}(J^{\dagger}J)_{n+1}\Delta]$ computing the acceptance would require just the action difference $\Delta S=\operatorname{Re}[S(x_{n+1}^{\prime})-S(x_{n}^{\prime})]$. This would require solving an implicit equation to determine the new configurations $x_{n+1}$. This is not the approach we will follow here. Instead, we follow a method based on the Grady algorithm Grady (1985); Creutz (1992) which can be designed to both be isotropic around $x_{n}^{\prime}$ and also to avoid explicit computation of $J$; we make proposals using

[TABLE]

where $\delta$ is a parameter used to adjust the step size so that the acceptance rate is reasonable. Note that this steps relies on the fact that the matrix $J^{\dagger}J$ is real, a property discussed earlier in the paper. The equation above involves $\left|\det\,J\right|$, but this will appear in our calculation only implicitly, since our algorithm will require only that we generate vectors $\eta_{\parallel}^{\prime}=J_{n}\Delta$ that are normally distributed; we only displayed the full probability function above to help us prove the detailed balance below. The acceptance rate is computed by first generating an auxiliary complex vector $\xi$ with probability

[TABLE]

and then the acceptance probability is computed using

[TABLE]

The total transition rate from $x_{n}$ to $x_{n+1}$ is then

[TABLE]

To prove the correctness of the method we follow a strategy outlined by Creutz Creutz (1992): we show that the product $P(x_{n})T(x_{n}\to x_{n+1})$ is symmetric in $x_{n}$ and $x_{n+1}$ which implies that the detailed balance is satisfied. We have

[TABLE]

This proves the correctness of our method. To derive the relation above we used the fact that $J^{\dagger}J$ is real and positive-definite, so that $\sqrt{\det J^{\dagger}J}=\left|\det J\right|$.

Our algorithm samples configurations $x\in{\mathds{R}}^{N}$ distributed with the desired probability $P(x)$ given in Eq. 20. The configurations are updated using the following steps:

Generate a complex vector $\eta^{\prime}\in{\mathds{C}}^{N}$ using the distribution $P(\eta^{\prime})\propto\exp(-\eta^{\prime\dagger}\eta^{\prime}/\delta^{2})$. 2. 2.

Compute $\Delta=\operatorname{Re}J_{n}^{-1}\eta^{\prime}$ and $\eta^{\prime}_{\parallel}=J_{n}\Delta$ which is the tangent component of $\eta^{\prime}$ at $x_{n}^{\prime}$. 3. 3.

Propose a new configuration $x_{n+1}=x_{n}+\Delta$. 4. 4.

Generate a complex vector $\zeta^{\prime}\in{\mathds{C}}^{N}$ using the probability distribution $P(\zeta^{\prime})\propto\exp(-\zeta^{\prime\dagger}\zeta^{\prime})$. 5. 5.

Compute $\xi=J_{n+1}^{-1}\zeta^{\prime}$ and use it to determine whether to accept the new configuration using ${{\rm P}_{\text{acc}}}$ in Eq. 23.

In step 1 and 4, $\eta^{\prime}$ and $\zeta^{\prime}$ can be generated by simply drawing each component from an appropriate normal distribution. In step 2, in order to show that $\eta^{\prime}_{\parallel}$ is tangent to $\mathcal{M}$ at $x^{\prime}$ we used the fact that we can decompose $\eta^{\prime}$ into components parallel and tangent: $\eta^{\prime}=\eta_{\parallel}^{\prime}+\eta_{\perp}^{\prime}=J(\eta_{\parallel}+i\eta_{\perp})$, where $\eta_{\parallel,\perp}\in{\mathds{R}}^{N}$, and thus $\Delta=\eta_{\parallel}=\operatorname{Re}J_{n}^{-1}\eta^{\prime}$, as it follows from the discussion surrounding Eq. 18.

The final ingredient necessary to implement steps 2 and 5 is a method to compute $J^{-1}v^{\prime}$ for an arbitrary complex vector $v^{\prime}\in{\mathds{C}}^{N}$ without first evaluating the matrix $J$. We note first that, for a real vector $v$, we can compute $Jv$ simply by integrating Eq. 13 to get ${\rm T}\!f(v)$, a calculation of complexity ${\cal O}(N)$. However, this does not work when the vector $v$ is complex, because the flow is non-linear when the vector has imaginary components, that is ${\rm T}\!f(v+iw)\neq{\rm T}\!f(v)+i{\rm T}\!f(w)$. For a complex vector the solution is to evolve the real and imaginary components separately, that it $Jv={\rm T}\!f(\operatorname{Re}v)+i{\rm T}\!f(\operatorname{Im}v)$. This requires two separate integrations of Eq. 13, but the complexity remains ${\cal O}(N)$. Armed with this routine, we can then compute $J^{-1}v^{\prime}$ for any vector by using an iterative method, for example GMRES Saad and Schultz (1986) or BiCGstab van der Vorst (1992). These algorithms allows us to compute $J^{-1}v^{\prime}$ through successive computations of $Jv$ without explicitly inverting $J$. We note that when the vector $v^{\prime}$ is tangent at $x^{\prime}$, its inverse $J^{-1}v^{\prime}\in{\mathds{R}}^{N}$ can be computed by integrating Eq. 13 backwards. While we did not use this property for this study, it is conceivable that this may be employed to optimize an inversion algorithm for general vectors.

The algorithm described in this section can be used to sample $P(x)$ without computing the matrix $J(x)$ and its determinant. To complete the calculation we need to evaluate the phase $\varphi(x)$. While there are methods that can estimate this phase accurately without computing $\det J$ Cristoforetti et al. (2014), in this paper we evaluate it directly by computing $\det J$ on the decorrelated configurations, which is a small subset of the one generated by the sampling algorithm.

Later on, in order to gauge the performance of our algorithm, we will compare it with a variant of our algorithm. For that purpose we improve the way proposals are made in relation to the method used in Alexandru et al. (2016b). The proposals are made according to the distribution ${\rm Pr}(x_{n}\to x_{n+1})\propto\exp(-\Delta^{T}J^{\dagger}(0)J(0)\Delta/\delta^{2})$ and accepted with probability ${{\rm P}_{\text{acc}}}=\min\{1,\exp(-\Delta S)\}$. This in effect samples the configurations not with the desired probability but with probability $P_{0}(x)\propto\left|\det J(0)\right|\exp(-\operatorname{Re}S(f(x)))$ so, when computing he observables the difference between the desired probability and $P_{0}(x)$ has to be reweighted by using

[TABLE]

where the average $\left\langle\cdot\right\rangle_{P_{0}}$ is taken with respect to the weight $P_{0}(x)$. For the model discussed in this paper, at weak coupling the matrix $J(x)$ will be close to $J(0)$ and the method should work well. As we increase the coupling the reweighting factor $\tilde{\Phi}$ oscillates quickly and the statistical power of the ensemble generated by $P_{0}$ will decrease making the process inefficient. We note that for this method we need to compute $J(0)^{-1}v$ for a large numbers of vectors. For sufficiently small systems, as the ones we consider in this paper, this can be done efficiently by precomputing $J(0)^{-1}$, allowing us to generate a large number of configurations very cheaply.

V Results

The data presented in this section refers to simulations that are run on lattices with $N_{t}=8$, $N_{\beta}=2$, and $N_{x}=8$. The total number of degrees of freedom for this system is $(2\times N_{t}+N_{\beta})\times N_{x}=144$. We set $m=1$ and the lattice spacing $a=0.2$. We run simulations for three different values of $\lambda$: $0.1$, $0.5$, and $1.0$. For each $\lambda$ we did short runs to determine the average sign $\operatorname{Re}\left\langle\Phi\right\rangle_{P}$ and increased $T_{\text{flow}}$ until it reached a value of around $0.2$. The parameters for these simulations, including the flow time $T_{\text{flow}}$, the average sign, and the number of updates for both our current algorithm and the old one are collected in Table 1.

The first issue we will address is the performance of the estimator that appears in the acceptance rate in Eq. 23. We can find the acceptance rate by integrating Eq. 22 over $\xi$; the result is a ratio of determinants $\det(J^{\dagger}J)_{n+1}/\det(J^{\dagger}J)_{n}$. One potential problem would be a low acceptance rate or, equivalently, the need to use very small proposals $\delta$. In Fig. 2 we show the relative size of the three components that enter in the acceptance rate: the action change $\Delta S$, the $\xi$-estimator $\xi^{\dagger}[(J^{\dagger}J)_{n+1}-(J^{\dagger}J)_{n}]\xi$, and the measure change $\Delta^{T}[(J^{\dagger}J)_{n+1}-(J^{\dagger}J)_{n}]\Delta$. We see that the action change dominates both the $\xi$-estimator and the measure change, for all values of the coupling. We conclude that there is very little loss associated with using the $\xi$-estimator rather than the determinant ratio in the acceptance rate.

The observable that we will focus on in this section is the correlator

[TABLE]

Note that the field $\phi$ is real thus $\phi(p)^{\dagger}=\phi(-p)$. We use periodic fields in the spatial direction so the momentum is quantized in units of $2\pi/N_{x}$. To compute the correlator $C(t,p)$ we compute the averages on the forward time leg of the Schwinger-Keldysh contour. To boost our statistics, we can also use the fields on the backward leg of the contour. Note that the correlator

[TABLE]

can be evaluated using the fields on the backward leg due to the reverse time ordering of the operator product. We assume here that $t\geq 0$. This correlator is related to $C$ via the relation $C(t,p)=\overline{C_{b}(t,-p)}$. Thus the observable we use is the following

[TABLE]

where $T_{b}=N_{t}+1+N_{\tau}/2$ is the beginning of the backward leg in our discretization. The formula above exploits the time translation invariance along the forward/backward legs. Note that this symmetry is exact only in the continuum: our discretization breaks this slightly due to the corner effects.

In the left panel of Fig. 3 we compare the results of our simulations with the perturbative calculation in the weak coupling ($\lambda=0.1$) and $p=0$ case. We see the the results from our two algorithms agree with each other—the error bars are smaller for $J_{0}$ because the number of statistics is significantly larger—and they agree very well with the perturbative results. Note that for this coupling the perturbation theory seems to work well: the effect of higher order terms gets smaller. For larger values of $\lambda$, however, the perturbative results should become unreliable. We quantify the convergence of perturbation theory on the right panel of Fig. 3 where we plot the value of $C(t=0,p=0)$ as a function of $\lambda$ at orders $\lambda^{0},\lambda^{1}$ and $\lambda^{2}$. We see that the series becomes unreliable well before we get to $\lambda=0.5$, with the second order term having larger magnitude than the first order correction, and the series diverging sharply away from our simulation points at $\lambda=0.5$ and $1.0$. We conclude that $\lambda=0.5$ and $1.0$ simulations are in the strong coupling region.

In Fig. 4 we plot the results for the real part of the correlator for all values of $\lambda$ used in this study. Since the statistics for the simulations using the Grady algorithm are smaller, the error bars are larger, so we only display the results for the largest coupling $\lambda=1.0$ (the ones for $\lambda=0.1$ are included in Fig 3.) We see that the correlator is systematically drifting away from its $\lambda=0$ value as we increase the value of the coupling. We have run statistical tests and determined that the results of the Grady and $J_{0}$ algorithm are statistically compatible. The apparent discrepancy in Fig. 4 is not statistically significant since the data-points for different time separations are strongly correlated.

In Fig. 5 we show the correlator for the smallest non-zero momentum, $p=2\pi/N_{x}$. We see that these correlators tend to have smaller statistical errors (this is true for higher momenta too), and they change only slightly as we increase the strength of the coupling. In the right panel of Fig. 5 we show the expectations from perturbation theory for this correlator: we see that the first and second order effects are much smaller than for the zero-momentum case, in agreement with the results of our simulations.

Before we conclude, we note that the performance of the “$J_{0}$” algorithm has been better than we anticipated. This was mainly due to the fact that for the parameters studied in this paper the fluctuations of the Jacobian were small. Thus the reweighting from the probability distribution $P_{0}$ sampled by the $J_{0}$ algorithm to the distribution $P(x)$ sampled by the Grady algorithm was very successful. To quantify this we measure the statistical power for each of the ensembles generated with $J_{0}$ algorithm

[TABLE]

where $N_{\text{cfg}}$ is the number of configurations in the ensemble, and $\phi_{i}$ are the configurations. This quantity is equal to one when all configurations contribute equally, which happens when $|\det J(\phi_{i})|$ does not fluctuate at all, and in the worst case it is $1/N_{\text{cfg}}$ when one configurations has a dominant contribution to the reweigthed ensemble. The statistical power of each ensemble is listed in Table 1. We see that as expected the statistical power decreases as we increase $\lambda$, but the its value is still $0.68$ even on the ensemble with $\lambda=1.0$, so that we can easily reweight.

VI Discussion and conclusions

We have explored two different algorithms to compute the path integrals arising in the Schwinger-Keldysh formalism. They are both based on deforming the contour of integration from real variables to a submanifold of the complexified field space. They both bypass the most difficult and costly part of the computation, namely, the calculation of the jacobian of the parametrization of the deformed submanifold by real variables. They also lead to more efficient, isotropic Monte Carlo proposals.

The first algorithm has general applicability and can be seen as an adaptation of the Grady algorithm previously proposed to deal with the fermion determinant. The second (“$J_{0}$”), uses a free field approximation of the jacobian (with the difference between the correct and the free field jacobian reweighted during measurements). This last algorithm breaks down at strong enough coupling but we observed that it performs efficiently well before the point where perturbation theory is no longer valid.

The algorithms were applied to the computation of real time thermal correlators in the $1+1$ dimensional $\phi^{4}$ scalar theory. The two methods agreed with each other and with perturbation theory results at small enough values of the coupling. The “$J_{0}$” algorithm is very efficient and its success at even relatively large values of the coupling is somewhat surprising. These calculations are, to our knowledge, the first reliable Monte Carlo real time calculations in a field theory.

The algorithms developed in this paper paves the way for larger scale calculations with finer lattices and/or larger number of spatial dimensions. The extension of the maximum time (bound in the present paper by $4\beta$) is a little more subtle and we uncover some evidence of trapping of the Monte Carlo chain in local minima of the effective action. In this case, the methods advocated in Fukuma and Umeda (2017); Alexandru et al. (2017b) should be useful and should be incorporated.

Acknowledgements.

A.A. is supported in part by the National Science Foundation CAREER grant PHY-1151648 and by U.S. Department of Energy grant DE-FG02-95ER40907. A.A. gratefully acknowledges the hospitality of the Physics Departments at the Universities of Maryland and Kentucky, and the Albert Einstein Center at the University of Bern where part of this work was carried out. P.F.B. and G.R. are supported by U.S. Department of Energy under Contract No. DE-FG02-93ER-40762. G.B. is supported by U.S. Department of Energy under Contract No. DE-FG02-01ER41195.