Real-time quantum dynamics, path integrals and the method of thimbles

Direct numerical evaluation of the real-time path integral has a well-known sign problem that makes convergence exponentially slow. One promising remedy is to use Picard-Lefschetz theory to flow the domain of the field variables into the complex plane, where the integral is better behaved. By Cauchy's theorem, the final value of the path integral is unchanged. Previous analyses have considered the case of real scalar fields in thermal equilibrium, employing a closed Schwinger-Keldysh time contour, allowing the evaluation of the full quantum correlation functions. Here we extend the analysis by not requiring a closed time path, instead allowing for an initial density matrix for out-of-equilibrium initial value problems. We are able to explicitly implement Gaussian initial conditions, and by separating the initial time and the later times into a two-step Monte-Carlo sampling, we are able to avoid the phenomenon of multiple thimbles. In fact, there exists one and only one thimble for each sample member of the initial density matrix. We demonstrate the approach through explicitly computing the real-time propagator for an interacting scalar in 0+1 dimensions, and find very good convergence allowing for comparison with perturbation theory and the classical-statistical approximation to real-time dynamics.

A Loop corrections 29 1 Introduction For a quantum system evolving from one state to another, Feynman's path integral quantization asserts that all possible paths in field space contribute to the quantum amplitude.These contributions are equal in magnitudes but have different phases [1].This poses a great challenge when one wants to compute the path integral through numerical methods, as although the interesting physics is often concentrated in some region of space of paths, a detailed cancellation of quickly oscillating functions must be achieved.The challenge is also known as the "sign problem", and appears whenever the path integral kernel cannot be made real by Wick rotation to a Euclidean action, such as when a chemical potential is present, or the correlators one is trying to compute involve a real time separation.
Recently in [2,3], it was shown that the real-time path integral can be computed through a Generalized Thimble Method, based on complexifying the field variables.In the complexified field space, one can deform the integration cycle of the path integral into the complex plane and still obtain the same result of the integral, so long as the integrand is holomorphic in the new complex variables.There is some freedom in how one deforms the contours, but a natural choice is to use a gradient flow (to be described below), starting from the original real field space.If we stop the flow at some finite flow-time, then the original field space will have flowed to some new field space M, and the integral over M is equivalent to the integral over original real field space.In particular, when the flow time approaches infinity, one finds that M ∞ is composed of Lefschetz thimbles [4], where the phase of the integrand is constant and therefore the "sign problem" is eliminated on each thimble (except for milder contributions from the residual phase [5]).In practise, we are not able to perform infinitely long flows.However, as long as the flow-time is large enough then the "sign problem" will be alleviated, in the sense that the highly oscillatory integrals of a function with constant magnitude (i.e. e iS/ ) turn into integrals of an oscillating function with decaying amplitude.
The Lefschetz thimbles are a set of special submanifolds within the complexifed space that contain critical points of the action, and points on the thimble will flow to (or from) these critical points.Equivalently, Lefschetz thimbles are the manifolds that are generated by the gradient flow from critical points.When there are many critical points it can be difficult, in practise, to determine which thimbles should contribute to the integral.The Generalized Thimble Method takes a finite flow time from the original integration manifold to M, and this manifold will approach the appropriate set of thimbles for the integral as the flow time is increased.Although this automatically selects the correct thimbles, in practise the numerical sampling algorithm can get stuck on one particular thimble, as the connections between the thimbles are exponentially small.This manifests itself as a multimodal problem in the Monte Carlo calculation of the integral when there is more than one thimble.Ideally then, one would prefer to work with systems that have a single thimble, and so a single critical point of the action.
We are interested in applying the thimble approach to real-time quantum systems, and in this context the critical points correspond to classical trajectories that extremize the action.The idea we shall follow, that allows us to work with a single thimble at a time, is that there is a single classical solution for a given initial position and velocity or, in the language of fields, a given ϕ(t = 0, x) and φ(t = 0, x).Of course, given that we are studying a quantum system, there will be an ensemble of initial positions and velocities described by an initial density matrix, but we will see that we are able to separate the path integral into a two-step sampling procedure; for each member of the initial condition ensemble, we may compute a well-defined contribution to the path integral using the Generalized Thimble Method, and subsequently average over the initial condition ensemble in a straightforward way.
The framework where one can separate the full path integral into an initial distribution and the subsequent dynamical part of the path integral already exists, and is known as the Schwinger-Keldysh, or in-in formalism [6,7].It is adapted to situations where one has initial data, rather than comparing in and out states, and as such one uses a time contour that starts at t 0 , extends to some T , and then goes back to t 0 , rather than to infinity.The value of T is arbitrary, so long as the path encompasses any operators O(t) that one is interested in.For some theoretical situations it is useful to take T → ∞, but for numerical simulations, such as in this paper, we work with finite T .We shall show in section 3.5 that it is the same reasoning behind the freedom of choosing T that enables us to separate the full path integral into two steps.
We will see that our approach to solving for the complete real-time quantum dynamics may be linked to popular approximation schemes, such as the classical-statistical approximation, (truncations of) real-time Schwinger-Dyson (Kadanoff-Baym) equations [9] or a quantum "dressing" of the classical path by Langevin methods in stochastic quantization [10].As for traditional Euclidean equilibrium lattice simulations we may compute the path integrals exactly from first principles, up to lattice discretization errors and finite numerical resources.
The structure of the paper is as follows: in section 2 we describe the Lefschetz Thimble Method and the Generalized Thimble Method, and introduce critical points and their role in evaluating the path integral.We connect to earlier work [3], and recall the flow equations to be used later on.In section 3 we discretize the path integral and show how the initial conditions may be separated from the remaining degrees of freedom, and set up the two-step sampling procedure.We set up a convenient parametrization of the discretized path integral variables entering in the real-time, but not necessarily closed-time path, path integral.We demonstrate how splitting up the sum over paths into subsets with fixed initial conditions, can resolve the multimodal problem in a straightforward way.We then explicitly derive the Gaussian initial density matrix, at finite temperature and in the vacuum, and take care of some technical points that arise.In section 4 we present our numerical model and algorithm for a field theory in any dimension, and demonstrate our approach for a theory in 0+1 dimensions, so quantum mechanics.We conclude in section 5. Some details of the perturbative one-loop correlator are placed in appendix A.

The path integral deformed into the complex plane
Consider the path integral written in the form1 , with real variables ϕ i , and I is a function of all ϕ i .Here we combine space-time indices into i, and will specify them more precisely later.As in the Feynman path integral, the exponent could be purely imaginary, so that the integrand is oscillatory with a constant amplitude.We can improve the convergence of the integral through complexifying ϕ i and, because of Cauchy's theorem, we can deform the real integration cycle into the complex plane and still obtain the same result for the integral.In the following we shall use ϕ i to denote the real field, and φ i shall denote the complexified field.As such, the initial integration manifold is R n , parametrized by ϕ i .This integration cycle is then deformed to a surface in C n with n real dimensions, parametrized by φ i .

Lefschetz Thimble Method
Such an approach is pioneered in [4,14] for Feynman's path integral, with the altered integration cycle known as Lefschetz thimbles, obtained by gradient flow, from critical points that are determined by ∂I/∂φ i | crit = 0.The over-line above refers to complex conjugation, and I is now considered a holomorphic function of the complex φ i .The Lefschetz thimbles are n-dimensional integration cycles in the n-dimensional complex (so 2n real dimensional) plane.As we can see from the flow equation, dI/dτ = i |∂I/∂φ i | 2 , Im[I] is constant on each thimble, and of the same value as at the critical point, while Re[I] keeps increasing with τ as we move away from the critical point, so its contribution to the integral (2.1) is exponentially suppressed away from the critical point.As a result, we achieve quick convergence by performing the integral on the Lefschetz thimbles.
The idea of integrating over Lefschetz thimbles can be naturally adopted to numerical simulations [5], especially through Monte Carlo methods with, for example, Langevin dynamics [15,17] and also Metropolis algorithms [16,21,22].In the following sections, we will use the term Lefschetz Thimble Method to refer to the methods of generating samples on Lefschetz thimbles.In the case of a single integration variable, the constraint that Im[I] is the same as it is at the critical point can almost determine the thimbles entirely [17].With more integration variables, however, this one constraint is not sufficient and we should return to using the gradient flow (2.2).To be precise, we should consider the flow starting from a small neighbourhood of the critical point on each contributing thimble, as the gradient flow will actually take infinite time to run away from the critical point itself.The neighbourhood should also be small enough to use an expansion of I up to quadratic terms, and with only these quadratic terms present we can solve the flow equation explicitly.This requires each isolated critical point, p, to be non-degenerate [5], ∂I ∂φ i p = 0, and det By Morse theory/Picard-Lefschetz theory, the matrix of second order derivatives of Re[I] has n positive eigenvalues and n negative ones, and near the critical point, we can approximate the Lefschetz thimble with the manifold generated by these n positive eigenvalues/eigenvectors.The "sign problem" is milder on the Lefschetz thimbles than on the real space.On each thimble, Im[I] is constant, and the only varying complex phase comes from the Jacobian of the transformation that maps the complex integration variables into real ones [16].More importantly, the exponential suppression of the magnitude away from the critical point makes the Monte Carlo simulation on each thimble possible.
A subtlety arises when there exists multiple critical points since in this case one has to find all the critical points and related Hessian matrices analytically, and then decide which combination of thimbles is equivalent to the original integration contour.There could exist one or more dominant critical points, giving similar contributions to the path integral.But there also exist concrete examples where equally dominant critical points cancel each other out in the integral, so that the main contribution comes from sub-dominant critical points [18,19].One might also have to sum over all contributing thimbles to not miss something [20].This is not an easy task for a general theory, so is there a technique that includes the complete integration cycle automatically, without having a "sign problem" at the same time?Such a technique is the Generalized Thimble Method.

Generalized Thimble Method
The gradient flow (2.2) serves two purposes.On the one hand, starting near critical points, it defines the corresponding Lefschetz thimbles.On the other hand, it maps the real integration cycle to the combination of thimbles contributing to the original integral.For instance, at τ = 0, we have the original n-dimensional real space, and as τ → +∞ we obtain the right combination of Lefschetz thimbles.In fact, the flow equation (2.2) in this case generates a family of n-manifolds that are characterized by the flow time, τ , and at any such flow time the integral would return the same result.Given the "sign problem" at τ = 0 and its absence at τ = +∞, one might expect the "sign problem" to be alleviated gradually along τ , and even at some finite τ the Monte Carlo simulation may already become effective.This turns out to be the case and such a finite τ approach, which is known as Generalized Thimble Method [24][25][26], has many applications in dealing with the "sign problem" in different scenarios [2,3,[23][24][25][26].
The finite τ manifold, M, has n real dimensions and is embedded in an n-dimensional complex plane.We can parametrize it with real variables as follows.Provided with initial real values ϕ i , the flow equation (2.2) transforms the fields into complex φ i .Thus we arrive at the equalities, The first equality is where we complexify ϕ i → φ i and perform the integration over the manifold M2 ; the second equality is where we think of φ i (τ f inal ) as a function of the initial ϕ i = φ i (τ = 0), and perform a co-ordinate transformation back to ϕ i .Note that in the final expression, I is evaluated at φ i (ϕ), while the first expression is evaluated at ϕ i , with the Jacobian providing the appropriate correction factor.
We can also deduce from the flow equation (2.2) that the Jacobian matrix J ij = ∂φ i /∂ϕ j satisfies d dτ with J ij an n × n identity matrix at τ = 0.In practice, one can carry out importance sampling with the weight P (ϕ) = e −Re[I]+ln |det(J)| , and then reweight by the remaining imaginary parts, We see this by noting that expectation values for operators are given by the following path integral While the Lefschetz Thimble Method approach is well-suited to an analytic approach, the Generalized Thimble Method with finite τ is more numerically oriented.On the other hand, the Generalized Thimble Method is not sensitive to the degeneracy of critical points.
To alleviate the "sign problem" one may have to go to a manifold with large τ , where the connection among different regions of the integration contour, flowing from multiple critical points, becomes exponentially small.As a result, simple Monte-Carlo sampling algorithms may get stuck in one region.This "multimodal" problem is a likely feature of the Generalized Thimble Method.Many sophisticated methods have been proposed to get the correct exploration of the manifold [23,27].But there is no doubt that both the Lefschetz Thimble Method and the Generalized Thimble Method are effective in the case of a single critical point.Then a natural question is whether we can tell the number of critical points beforehand.It turns out that we can, at least for a scalar theory.

Theoretical developments for the real-time path integral
At this point, we will derive a series of results for the path integral, which will all come into play, when we put together our algorithm in section 4.

The path integral
To fix our conventions we will start by deriving the path integral expression for calculating operator expectation values in the Heisenberg picture, Ô Φ, Π , with operator Ô consisting of the scalar field operator Φ and its canonical conjugate, Π, at one or more times.
We follow the convention of [11] : where |φ; t and |π; t are eigenvectors of operator Φ(t) and Π(t) respectively.In the formulae above, a discretized d-dimensional space was assumed.That is, N s sites along each spatial direction and distance dx between two neighbouring sites, so the volume V = (N s dx) d and, furthermore, we suppressed the spatial index.For instance, Dφ = x dφ(x).
It is also convenient to switch between continuous and discrete expressions via, We can then calculate Ô Φ, Π by inserting complete sets of |φ; t i φ; t i | in succession along the temporal direction, leading to with ρ Φ(t 0 ), Π(t 0 ) the initial density matrix operator at t 0 .Fig. 1 gives a graphic demonstration of the insertion along the temporal direction.
Here we separate φ + and φ − vertically for demonstration purpose.All these fields live on the real-time line.The difference between two neighbouring t is a constant, given by dt.
In the presence of operators O, the insertion is not unique.There are two features worth noting.(1) There are different ways for the operators to appear in the expression.For instance, in the case of Ô = Φ(t α ) Φ(t β ) and t α > t β , if Φ(t β ) appears in the upper (φ + ) layer, then Φ(t α ) can appear either in the upper (φ + ) or lower (φ − ) layer.We will see what this implies for the path integral in section 3.9.(2) One is free to choose the turning point φ m , as long as the contour includes the operator Ô.The path integrals with different turning points give the same expectation value of the operator.
First, we need to calculate each Feynman kernel φ i ; t i |φ j ; t j .Here we only assume that the time difference |t i − t j | is small, but do not specify which time is earlier.Since we want to derive the path integral with dt finite, a symmetric expression of the kernel seems a better choice, as it will converge more quickly in the limit dt → 0. Thus by evolving each state to the equal time t = (t i + t j )/2, and then inserting the complete set of |π; t π; t|, we arrive at the expression, where the operator Ĥ is the Hamiltonian, which contains only up to quadratic terms of Π.
For the scalar theory, we assume the general expression, with Ĉ( Φ) composed of spatial derivative terms and a field potential.We do not need to know the exact expression of Ĉ( Φ) at the moment, but demand Ĉ( Φ) is local in time.We also assume that all operators, for instance Ĥ, may be written as functions of variables φ and π.The function H(φ i , π) is then the result of the operator Ĥ acting on states.Given the Hamiltonian, the Lagrangian is which is symmetric in i and j.In light of eqn.(3.4), the wave function φ; t|in = Dφ φ; t|φ ; t − dt φ ; t − dt|in satisfies the Schrödinger functional equation [13], in the limit dt → 0. Thus, Feynman's kernel is the propagator for small time intervals.We emphasize that the derivations in (3.4) are valid for both t i > t j and t i < t j .From φ + m−1 to φ m , the time difference is dt, but from φ m to φ − m−1 , it is −dt.Now we can continue working with eqn.(3.3) where N is a collection of numerical constants that appear in kernel (3.4), and the integration contour C is understood as the contour shown in fig. 1.In the discrete theory, the integral over C in the exponent is really an abbreviation of, where, to write the expression elegantly, we denote φ m = φ + m = φ − m .On the other hand, since the numerical constant N does not depend on the operator Ô, we can fix it by taking the case Ô = 1, where we utilize the fact that the trace of the density matrix is one.Therefore, we can write the expectation value of the operator as,

Critical points
We are now in a position to find the critical points in eqn.(3.12).We write I = −i C dtL/ + • • • , with ellipsis denoting extra terms coming from the initial density matrix, which are only functions of φ + 0 and φ − 0 .To study the critical points it is convenient to use another basis, φ cl and φ q , defined through [6-8, 29, 30] 3 , With these 4 , the action becomes 3 In the literature, there exist alternative ways to transform φ + and φ − , with Keldysh's original convention [7,30] corresponding to φ ± = φ cl ± φ q / √ 2.Here we follow the approach of [8,29], but we adopt the names φ cl and φ q from [30]. 4 Even though we do not apply the change of basis to φm(x), as there is only one field, it will be useful to introduce φ cl m (x) = φm(x) and φ q m (x) = 0.But we do not treat φ q m (x) as a variable. where We may derive two general results without knowing the explicit form of the Lagrangian: 1. E m = 0.The only term in the exponent containing φ m (x) is the product of φ m (x) and φ q m−1 (x).Actually, in eqn.(3.12), one can integrate φ m (x) out, and get a delta function, as follows, If one further integrates out φ q m−1 (x), eqn.(3.12) would become the same form as the original integral, but with the turning point φ m replaced by φ cl m−1 , and with an extra overall constant.We emphasize the fact that the integration over (φ cl m , φ q m−1 ) together is a constant, and it will not alter the remaining path integral, except through the overall constant.One may integrate out the (φ cl i , φ q i−1 ) one by one, as they become the last pair along the real-time direction.By continuing this process down to φ 0 , we arrive at This is just eqn.(3.10), written in reverse order, and also provides an alternative way to compute the constant N .Of course, to avoid keeping numerical constants, one can execute such contraction simultaneously in both the numerator and denominator of eqn.(3.11).However, the contraction in the numerator is no longer valid once Ô(t) is reached.Generally, if t max is the maximum time that the operator Ô depends on, then as long as t m t max , the path beyond t max is contractible.This corresponds to the freedom that one can have in choosing the closed time path when restricted to the real-time line.
2. All terms in E i contain odd powers of φ q i (x), as even powers of φ q i (x) cancel out.One can check this by expanding eqns.(3.14) and (3.15) as a Taylor series in φ q i .In fact, the quantum field theory can be computed in perturbation theory of φ q [8].The leading order theory has a term linear in φ q appearing in the exponent, and if we carry out the integration of φ q explicitly, the leading order theory is the classical theory.A simple example is λφ 4 theory (suppressing for moment the initial density matrix part of the expression), By keeping the leading term in the final factor, and then integrating out φ q , we find the delta function, which means, in the leading order theory, that φ cl satisfies the equation of motion of the classical field.More generally, ∂I ∂φ q φ q =0 = 0 leads to the classical equation of motion.Furthermore, when φ q i (x) = 0 at any x, then ∂E i /∂φ cl i (x) must also vanish, since it consists of odd terms of φ q i .
We may write down straightforwardly for 0 and We now note that the critical points are determined by ∂I/∂φ| crit = 0 for all φ, from which it follows that eqns.(3.19)-(3.21)all vanish at those points.We can now show by induction, that critical points require all φ q i (x) = 0 with 0 < i < m.This is true for i = m − 1, as the vanishing eqn.(3.21) alone indicates φ q m−1 (x) = 0 at any x.Furthermore, if φ q i+1 (x) = 0 along with φ q i (x) = 0 at any x, then as this implies ∂E i /∂φ cl i (x) = 0, we see that the vanishing of eqn.(3.20) leads to φ q i−1 (x) = 0. We can apply this induction down to ∂I/∂φ cl 2 = 0, such that all φ q i (x) = 0 with 0 < i < m.Now that we have φ q i (x) = 0 at the critical point, we can use the vanishing of eqn.(3.19), i.e. ∂I/∂φ q i (x) = 0, to lead us to the classical equation of motion, Notice that the second term on the left-hand side contains only φ cl i .Therefore, eqn.(3.22) determines φ cl i+1 (x) uniquely once φ cl i (x) and φ cl i−1 (x) are known.In other words, once φ cl 0 (x) and φ cl 1 (x) are known, we can uniquely solve all subsequent φ cl .In this sense, we can assert that the critical points are completely determined by φ cl 0 (x) and φ cl 1 (x).What we have shown, therefore, is that there is a single critical point for each given φ cl 0 (x) and φ cl 1 (x), and so by picking φ cl 0 (x) and φ cl 1 (x) there will be a single thimble associated to that single critical point.We now need a scheme to select φ cl 0 (x) and φ cl 1 (x), and for this we need an explicit expression of the initial density matrix.

Thermal initial density matrix for a free field
Later on, we will be particularly interested in Gaussian initial conditions, which may then be chosen to be vacuum, thermal equilibrium or any out-of-equilibrium initial Gaussian state.
But before we specialise to Gaussian states, we will first recall how a general thermal equilibrium state may be introduced as a path integral of imaginary time.
The density matrix operator for thermal equilibrium is ρ = e −β Ĥ /Z, where 1/β = k B T , with k B being Boltzmann's constant and T the temperature.The normalization Z = Tr e −β Ĥ is just an overall constant, which we will suppress for now.In this case, the insertion of complete sets leads to, with dβ = β/N .As the label suggests, it would be convenient to also denote φ + 0 as φ 0 and φ − 0 as φ N .The computation of each single kernel is similar to (3.4), and we can also compute it in a symmetric way, where the Lagrangian is defined similarly to the real-time one, but with dt substituted by It is then straightforward to compose the expectation value as a series of integrals, along a trajectory from φ N (so φ − 0 ) to φ 0 (so φ + 0 ), through negative imaginary time, In combination with the integral along the real-time as in eqn.(3.12), the whole path integral is defined on a closed contour in the complex time plane, which is periodic along the imaginary time, with a period β.Since there exist different ways to insert complete sets, there is some freedom in choosing the contour in the complex time plane.For a graphic illustration, see fig. 2. So far, we have considered the density matrix of a general scalar field, but for free fields we can carry out the integrals in eqn.(3.26).It is more convenient to do this in momentum space, so that we introduce For thermal equilibrium, the complex time path is periodic along the imaginary time direction, with the period β, and there is some freedom in choosing the contour in the complex time plane.(L) The Schwinger-Keldysh closed time contour used in [2,3]; (R) The Schwinger-Keldysh closed time contour used in [8].In section 3.3 we use the right-hand side path to derive analytic expressions, where both trajectories of t 0 → t m and t m → t 0 are located on the real-time line, and the vertical offset between them exists only for demonstration purpose.
Since φ is a real field, φ(−p) = φ(p) † , and we may write Thus it would be more appropriate to use its real and imaginary components as integration variables, in particular √ 2φ re (p) and √ 2φ im (p), which can be regarded as the result of a unitary transformation of (φ(p), φ(−p)).On the other hand, one can also arrive at the same variables, by performing a real-to-real Fourier transform in the first place.Later on, we will use p, re, im to mean that it is these real integration variables that we use.But it is easy to switch between (φ(p), φ(−p)) and ( √ 2φ re (p), √ 2φ im (p)).so that the free Lagrangian in momentum space takes the form, where ω p = p 2 + m 2 , and V is the spatial volume. 5We can now switch eqn.(3.26) into momentum space, and carry out the integrals, where the overall constant on the second line is changed due to the Fourier transform, and to reach the last line we take the limit dβ → 0. We are now able to calculate the partition function as, (3.31) 3.4 Initial density matrix for vacuum and n-particle states Alternatively, we can also derive everything from the n-particle eigenstates.The free theory is equivalent to a sum of independent harmonic oscillators with different ω p .Therefore, one can derive n-particle eigenstates for the free field theory as one does in the harmonic oscillator.We will skip the details of the derivation and only provide the final formulae.
In momentum space, the vacuum wave function is With it, we can write the density matrix of the vacuum state as, The wave function of the n-particle state is where the Hermite polynomial h n (z) is defined as: We can now compute the density matrix of any pure state or mixed state, as long as it can be expanded with n-particle states.For instance, it is straightforward to calculate the density matrix for the thermal states, up to the partition function Z, where to get the final expression, we have used Mehler's formula In the limit β → +∞, it becomes to (3.33).The density matrix of the thermal state at zero temperature gives the density matrix of the vacuum.So we are going to stick with the free thermal density matrix in the following sections, and treat the vacuum state as a special case.

Path integral with a free initial density matrix
Given a free initial density matrix, the full path integral has the general form, or, in the φ cl and φ q basis, with the occupation number given by The initial density matrix in (3.40) implies that the field φ cl 0 is drawn from a normal distribution with the variance proportional to 2n p + 1, while φ q 0 comes from a normal distribution with variance proportional to 1/(2n p + 1).We can get a better understanding of this observation by integrating out φ q 0 , noting that φ q 0 also appears in the last term of eqn.(3.40).However, by assuming the theory to be free at t 0 , we will not encounter any higher order terms of φ q 0 , i and we see that φ q 0 interacts only with φ cl 0 and φ cl 1 .After the integrating out φ q 0 the path integral takes the form, where L denotes L with all φ q 0 related terms removed.One now recognizes the new term in the square bracket above as just the time derivative of the scalar, but with finite dt, and we now see that the density matrix gives Gaussian distributions to φ cl 0 and φcl 0 with variances given by, In section 3.2, we mentioned that in the perturbation theory of φ q , the leading order theory has linear φ q terms in the exponent, and therefore one can integrate φ q out and obtain the classical equation of motion.There is still, however, the initial density matrix left.This means that the initialization of the classical theory should respect the distribution (3.45).
In practice, we can generate ensembles of initializations of φ cl 0 and φ cl 1 according to (3.44) and (3.45), and then use (3.22) to find the full classical history.As we will show below, this classical history may then be used as the starting point for our Monte Carlo simulation of the path integral, although the Monte Carlo process essentially washes out the memory of the classical history (except φ cl 0 and φ cl 1 , which are held fixed for a given Monte Carlo run.).
In the full quantum field theory, we also want to separate the initial density matrix contribution from the rest of the closed time path in the path integral.There are two reasons for doing this: (1) It is much easier to write the initial density matrix part in momentum space, and the subsequent dynamical part of the path integral in configuration space.
(2) There is no "sign problem" in the initial density matrix piece.In fact, the distributions in the initial density matrix piece of (3.45) are ordinary Gaussian distributions, and simple Monte Carlo methods are sufficient to generate samples of φ cl 0 and φ cl 1 .Thus, in addition, we also want to treat φ cl 0 and φ cl 1 on a different footing from the other integration variables.However, while the initial density matrix part involves only φ cl 0 and φ cl 1 , the remaining part of the path integral also contains φ cl 0 and φ cl 1 .So is the separation legitimate?The answer is yes, but with a note of caution.

Separating variables
When separating φ cl 0 and φ cl 1 from the other integration variables, we should check that the following equality is valid, where ρ φ cl 0 , φ cl 1 is the density matrix part in eqn.(3.43), and is a function of φ cl 0 and φ cl 1 only.Apparently, to have the equality valid, the lifted integral should be independent of φ cl 0 and φ cl 1 .To show that this is true, we make use of a feature that we have already explored: The only term in L containing φ cl m is from φ cl m (x)φ q m−1 (x), and by integrating out φ cl m , we obtain a delta function, δ(φ q m−1 ).Then by integrating out φ q m−1 , we obtain an integral similar to the previous one, but with φ cl m−1 now playing the role of φ cl m .We can continue this contraction of the closed time path down to φ q 1 , where we then find δ(φ q 1 ).Now, we know that all φ cl 0 and φ cl 1 appear in L only through their products with φ q 1 , so by integrating out the delta function of φ q 1 , we know the result has no dependence on φ cl 0 and φ cl 1 .Concretely, the result of the integral is which is independent of φ cl 0 and φ cl 1 , and so a constant from the point of view of the integral over initial conditions.We may thus perform the separation of variables in (3.46).

One critical point for one initialization
We separate the whole path integral into two parts: the initial density matrix and the rest of the path integral.To implement the Monte Carlo simulation, we propose different algorithms for each of these different parts.
1. We assume the initial density matrix is known, so we can sample φ cl 0 and φ cl 1 directly according to the initial density matrix, using simple Monte Carlo algorithms.There is no "sign problem" in the procedure, as in momentum space the distribution function is real and vanishes exponentially as |φ| → ∞ [12].Notice that the initial density matrix is a function of φ cl 0 and φ cl 1 only, but the rest of the path integral also depends on φ cl 0 and φ cl 1 .We denote such sampled fields as φcl 0 and φcl 1 , and a Fourier transform is necessary to bring the fields into configuration space for later use.All these φcl 0 (x) and φcl 1 (x) are real.2. Provided with each φcl 0 and φcl 1 , we then perform importance sampling according to with the Generalized Thimble Method, according to an algorithm such as in [3].Note that the quantum and classical fields start at 1 and 2 in the product, respectively, because φ q 0 has been integrated out, while φ cl 1 and φ cl 2 are specified as initial data for each initialization.The sampled φ cl i+1 and φ q i with 1 ≤ i ≤ m − 1 in this procedure are complex.With reweighting (2.6), we can calculate the expectation value of an operator Ô over a single initialization, which is equivalent to, The full expectation, Ô , in eqn.(3.46) will then be the mean of all the singles, Ô single .For the integral (3.48) above, we can repeat the analysis in section 3.2 to find all the critical points, this time with I = −i C dtL / .In fact, the conclusions in section 3.2 are still valid here: At critical points, all φ q i (x) = 0, so I = 0, as it consists of odd terms of φ q , and all φ cl i+1 (x) are uniquely determined through the classical equation of motion (3.22), once φcl 0 and φcl 1 are specified.In other words, for each initialization, there exists one and only one critical point.This means that for step 2 above, we will not encounter any multimodal problem that would be caused by the existence of multiple critical points.
However, the initial density matrix could possess multiple saddle points in its distribution.For instance, we expect this to happen in the density matrix of n-particle state when n = 0, or in the case of multi-scalar fields where there exists some symmetry among those scalars.Still, this will not change the conclusion that there exists one and only one critical point for the thimble part of the calculation, and we only need to deal with one thimble/critical point on step 2.
We stress that the derivation is valid on the complexified fields, and the thimble must contribute to the original integral, as the critical point is located on the real field plane.There is one more thing we can predict.With each initialization, the averaged phase e −iIm[I]+iarg(det(J)) P must be real and positive, due to eqn.(3.47).Furthermore, on the Lefschetz thimble, I vanishes at the critical point, so Im[I] = 0 on the whole thimble, and only the residual phase arg (det(J)) contributes.

Two-point functions
In order to test the formalism we will calculate the two-point correlators analytically, and compare them with numerical results based on the procedure described above.One can do this in the framework of perturbation theory, that is we first compute free correlators and then add the loop corrections.In this section, we only explicitly derive the free two-point functions, while a 1-loop correction will be included in App A. See also [8].Since in the free theory, different momentum modes are independent of each other, we can focus the calculation on a single mode.There are two equivalent ways, up to a constant due to the integration of φ q 0 , to write the path integral, with constants where ω p is the frequency in the continuous theory but, because of the discretization, it is ωp that propagates on the lattice.In the limit dt → 0, ωp converges to ω p .For finite dt, it is convenient to replace ω p in (3.50) and (3.51) with sin(ω p dt)/dt.With only Gaussian functions in (3.50) and (3.51), we can calculate the free two-point functions as, where A and x are understood to be a symmetric complex matrix and a real vector respectively.The size is given by the number of discrete points on the time contour of choice.The above normalization is appropriate for the discrete theory, while for the continuous theory, there will exist a factor of V in the definition.To compensate this, we simply assume V = 1 in the following derivation.

Time-ordered correlators
It is straightforward to identify the matrix A in eqn.(3.50), then calculate its inverse, and use (3.53) to discover that the two-point functions in the (φ + , φ − ) basis are where the star denotes complex conjugation, and the matrix F is There are two features worth emphasizing in the above expression.
1.In the vacuum, that is n p = 0, we notice that the rows and columns corresponding to φ + 0 → φ m (i.e. the upper-left part of F ) lead to F jk = exp (−iω p |t j − t k |), and give the Feynman propagator, which is defined as e −iωp|tx−ty|+ip(x−y) 2ω p .
(3.56) Thus we get the correct i prescription in the propagator.This also means the correlators φ + i φ + j 0 are time-ordered, while the correlators φ − i φ − j 0 are anti-time-ordered.On the other hand, when n p = 0, we can calculate the equal-time correlator through summing the Matsubara frequencies, This corresponds to calculating the equal-time elements in eqn.(3.54).
2. There exist symmetries in the above two-point functions.For instance, In fact, although we can have many integration variables φ i at time t i , there is only one operator Φi , and it is actually easier to discern the symmetries from the operator formalism, with On the other hand, as in eqn.(3.54), and this makes manifest the KMS condition G > (t i − t j ) = G < (t i − t j + i β) [8].

Classical-Classical and Quantum-Classical correlators
We could obtain the correlators such as φ cl i φ cl j or φ q i φ cl j through a rotation of φ ± i φ ± j in eqn.(3.54), but it is instructive to derive the expression from scratch with a simple example.Consider m = 3.Then the matrix A in eqn.(3.51) is with constants We treat φ m as a φ cl field.Since we have also integrated out φ q 0 , in the end there are two more φ cl fields than φ q fields.Following eqn.(3.53), we arrive at where This may be summarized by the following: φ q i φ q j 0 = 0, (3.68) We see, for example, that the correlators φ q i φ cl j vanish unless i < j, and so correspond to the advanced propagators.Furthermore, because of the advanced propagators, any loop correction will not alter φ q φ q = 0. Actually, we can derive this conclusion much more quickly from the operator formalism (3.58):

Numerical Simulation
We now demonstrate how to carry out numerical simulations, with an example of λφ 4 theory (see also [3,8]), using the following action, Ideally, we would like to simulate a 1 + 1 or even 3 + 1-dimensional system.But in those cases, one should stick with some specific renormalization scheme in order to compare with the result of continuum theory.This is beyond the scope of the present work, and is postponed for later work.Instead, we find it is straightforward to compare with theoretical predictions in 0+1-dimensional system, so quantum mechanics 8 , where no divergence exists, and therefore no renormalization scheme is required.We shall set up our definitions in d = 1 spatial dimensions, whereas in the actual simulations presented here, we have further reduced to d = 0 quantum mechanics.Throughout the paper, we set mdt = 0.75 for small couplings, and mdt = 0.5 for large couplings, (more details in our future publications).Space is discretized on N s sites, with periodic boundary conditions, and the time direction is discretized as above onto N t = 2m + 1 sites going back and forth on the Keldysh contour (see fig. 1).

Warm-up: Classical statistical approximation
We set the initial φ cl 0 (p) and φ cl 1 (p) according to eqn.(3.45), a Gaussian thermal density matrix. 9Given the distribution, we generate random samples of momentum-space variables φ cl 0 (p) and φ cl 1 (p) which are then Fourier transformed to position space φ cl 0 (x) and φ cl 1 (x).Now, we can compute the classical field evolution through the equation of motion, We use φ to refer to the fact that these are not variables of integration in the path integral.They represent the critical configuration in our complexified field configuration space, φ cl = φcl , φ q = 0, from which we will initiate our Monte-Carlo simulation in later sections.
Fig. 3 (left) shows the correlator for a single such classical trajectory.In a classical simulation, we can only compute the classical-classical correlator.By averaging over the ensemble of initial conditions, we recover the "classical-statistical" approximation to quantum dynamics, shown in fig. 3 (right).We show the results for a free field, λ = 0 and an interacting theory λ = 0.2.The correlators are very similar, but deviate enough that we can tell the difference with moderate statistics.The loop calculation is discussed in appendix A, where it is found that at 1-loop we just need to make the replacement ω 2 p → ω 2 p + λ 4ω .This is substituted into (3.52) to find ωp , which is then used in expression (3.65) for the classical-classical correlator.

Warm-up: Quantum average of a single initial realisation
Going beyond the classical approximation then amounts to performing the complete path integral, the integrations of all the field variables not associated with the initial conditon, It turns out that the exponent I is more conveniently expressed in the (φ + , φ − ) basis than using (φ cl , φ q ), as the interaction terms are simpler there.We therefore switch to (φ + , φ − ), except that at t 1 should be treated differently, since we count φ cl 1 into the initial condition, leaving φ q 1 as the only variable at t 1 .The exponent I also contains φcl 0 and φcl 1 , and may be written as where we have adopted a field redefinition as illustrated in fig.4, and the time differences are denoted as In the exponent, there are terms like φ q 1 (x) φcl 0 (x) − 2φ q 1 (x) φcl 1 (x) + • • • , where φcl 0 and φcl 1 can appear.In fact, an extra φcl 2 (x)φ q 1 (x) term will cancel out these linear-in-φ q 1 (x) terms, due to the equation of motion (4.2).Therefore, we are able to substitute these terms with φcl 2 (x) term only, and this simplifies expression (4.3) a lot.Given that φ cl 1 is part of the specified initial data, we define φ 1 = φ q 1 /2 to ensure that at site 1 only φ q 1 is included in the dynamical part of the path integral.To arrive at eqn. (4.3), we have also used that, There are N tot = N s (2m−2) variables in total, and we will adopt a more compact notation, merging space and time labels into a single integer a.For all the field variables φ a , we start our Monte-Carlo chain for the dynamical part of the path integral from φa , the classical critical-point configuration.In subsequent Monte-Carlo steps, these will be changed into new real values ϕ.For each such value, the gradient flow equation into the complex plane now reads The Jacobian matrix J itself, defined with element J ab = ∂φ a /∂ϕ b , evolves along the flow as, where a summation over index l is understood.
For the Lefschetz Thimble Method, then J(τ = 0) is determined by the eigenvectors of positive eigenvalues [5] of the Hessian evaluated on the critical point field configuration.We use the Generalized Thimble Method, then J(τ = 0) is just the identity matrix [3].With these flow equations, one can now apply thimble methods to generate samples for the dynamical part of the path integral.For more on algorithms based on the Lefschetz Thimble Method, see [5,15,16].And for more on algorithms based on the Generalized Thimble Method, see [21][22][23][24][25][26][27][28].
Our algorithm can be briefly summarized as follows: 1. Generate an initial value for φ cl 0 and φ cl 1 according to a Gaussian distribution given by eqn.(3.45).Determine the critical configuration by solving eqn.(4.2).
Evolve φ and J from τ = 0 to τ = τ f , for some final flow time τ f 3. To go from the n-th to the n + 1-th configuration in our Monte-Carlo chain for the dynamical part of the path integral, first propose the (n + 1)-th configuration ϕ n+1 = ϕ n + ∆, where the vector ∆ follows the proposal distribution, with some constant parameter δ10 .
4. Use the gradient flow equation to evolve φ n+1 and J n+1 from τ = 0 to τ = τ f . 5. Accept or reject new configuration according the acceptance probability If the new configuration is rejected, choose the (n + 1)-th configuration to be the same as the n-th configuration.
6. Repeat ( 3)-( 5) until we have enough statistically independent configurations to average over, for this one initial condition realisation.
7. Repeat ( 1)-( 6) for n initial times, to get enough initial conditions to average over (these are statistically independent by construction).
On the thimble approach (3)-( 5), we follow the prescription given by [3], with the difference that we perform an LU decomposition for matrix J to calculate its inverse and determinant directly.Therefore, we can have the acceptance probability with the explicit existence of det J.With the proposal distribution and acceptance probability above, the obtained samples will admit the probability weight P = e −Re[I]+ln |det J| .The numerical effort is substantial, and many technical details, performance tests and detailed numerical investigations will be reported in our future publications.
In fig. 5, we show the correlator for a single classical trajectory, and compare it to the correlator when averaging over the quantum variables (but without averaging over initial conditions, only step 1-6 of our algorithm).In the left-hand plot for the free theory (λ = 0), in the right-hand plot including interactions (λ = 0.2).We see that the quantum averaging is has only a small effect for the free theory, whereas including a moderate interaction strength there is statistically significant effect, increasing over time.

All warmed up: Full quantum evolution
We are now ready to carry out the inner (Monte-Carlo integration on the thimble) and outer (initial conditions) integration together, to find the full quantum correlator, given our initial Gaussian state.The simulations presented here use n initial = 200 ∼ 60 initializations, with (5 ∼ 20) × 10 5 Metropolis updates for single initialisation, in order to give small enough statistical errors.Fig. 6 (left) shows the two-point cl-cl correlator for the full classical-statistical simulation (pink) and the full quantum simulation (black).Overlaid also the 1-loop perturbative result (in red).Fig. 6 (right) arises from subtracting the free propagator, to highlight the contribution from interactions.We see that the classical-statistical approximation performs very well at these values of the coupling, and that apparently the differences arising from quantum averaging each initial condition (fig.5) are in turn largely washed out when averaging over initial conditions.The 1-loop approximation shown in red is distinct from the other two curves, showing that we are not in the extreme small-coupling limit, and so the agreement between classical-statistical and quantum approaches does apply to an interacting system.
We now proceed to increase the coupling λ, beyond the naively perturbative domain.We show in fig.7 the case λ = 4, where we can now clearly distinguish the classicalstatistical (pink) from the fully quantum result (black).They are both different from the free theory (green) and the 1-loop approximation (red).

Conclusions
Real-time quantum dynamics is well-defined in terms of the Schwinger-Keldysh formalism, and although the classical-statistical approximation often does very well in some cases, simulations of truncated Kadanoff-Baym equations have shown that quantum corrections are important in other contexts.
We have investigated a new, in principle exact, method for computing real-time quantum correlators directly from the path integral.This is possible through Monte-Carlo sampling, as the sign problem inherent to the complex action can be softened by flowing the field variables into the complex plane.We have presented a number of technical developments necessary to generalise the work of [2,3] to initial-value problems.For a discrete space-time, we have verified that the scalar field path integral can be separated into two parts: the initial density matrix and the following dynamical part.Under such a separation there exists one and only one critical point, which helps when we implement either the Lefschetz Thimble Method, or the Generalized Thimble Method on the dynamical part.We use a symmetric discretization of the theory, in both a symmetric Feynman kernel and a symmetric time contour.With such a discretization we can find all the critical points.
To demonstrate the implementation of our approach, we have computed the real-time propagator for a scalar field in 0+1 dimensions, with a Gaussian (free-field) initial condition.We found good statistical convergence, and agreement with the free analytic correlator (up to discretization errors).Once interactions were included and increased we found that we could distinguish from the free case, that the 1-loop perturbative result began to fail, and that for very large couplings, the classical-statistical approximation became unreliable.
In the present paper we have used the initial density matrix of the free theory, as in this case, we can integrate out φ q 0 explicitly, allowing us to obtain the familiar initial distribution of φ cl 0 and φcl 0 .There is no difficulty in extending the calculation to the case of a more general density matrix, as long as we know how to generate the initialization for φ cl 0 and φcl 0 .Note, however, that a density matrix containing φ q 0 and φ cl 1 might still be plagued with the "sign problem" owing to the appearance of a factor of iφ cl 1 φ q 0 in (3.42).This only affects the density matrix part of the path integral, so the thimble approach may still be used for the remaining dynamical part.On the other hand, we have also in mind that real physical situations can be modeled by turning on the interaction after the initialization, either instantly or gradually, and the method developed in the present paper can naturally deal with time dependent interaction coefficients.The computational cost of the thimble approach is aO(n 3 ), with n the number of variables and a the number of samples.By separating the simulation into two parts with n 1 and n 2 variables respectively, the cost becomes a 1 O(n 3 1 ) + a 1 a 2 O(n 3 2 ), corresponding to generating a 1 different initializations and for each initialization a 2 Monte Carlo samples.If a is not sensitive to n, the cost will be smaller than aO((n 1 + n 2 ) 3 ), when n 1 and n 2 are big numbers.In fact, if this is the case, we can further separate the path integral into more pieces, with each piece depending only on its predecessor but not successor, as each piece becomes an initial condition for the part that follows it.
We have postponed a number of numerical technicalities, diagnostics of the method and further numerical tests of various aspects of the approach to a future publication.Simulations on more general initial conditions and potentials, and in 1+1 dimensions are also underway.
Acknowledgements: PMS would like to thank Alexandru for some email correspondence.ZGM would like to thank Prof. Bedaque for useful suggestions on the algorithm.PMS and SW were supported by STFC Grant No. ST/L000393/1 and ST/P000703/1.AT and ZGM are supported by a UiS-ToppForsk grant.The authors were also supported by a ECIU travel grant.The numerical work was performed on the Abel supercomputing cluster of the Norwegian computing network Notur.

A Loop corrections
In this section we shall look at the loop corrections to the two-point functions, and we shall be using the continuum expressions in order to provide approximate expressions to the discrete case.First we look at the loop corrections to the Feynman propagator, and then we will see how the computation is adapted to the (φ cl , φ q ) basis.
The Feynman propagator is given in (3.56) as i dω 2π e −iω(tx−ty )+ip(x−y) ω 2 −ω 2 p +i , while the interaction vertex is − iλ 4! .The loop correction to the propagator is shown in Fig. 8, where the thick solid lines correspond to the Feynman propagator.This may be calculated in zero + +… spatial dimensions as follows.
where δm 2 = λ 4ωp , and we have used It is also instructive to use the (φ cl , φ q ) basis, for which we shall use the continuum expressions to give an approximation to the discrete calculation, and so we start by noting from (3.65-3.68)that the continuum propagators are given by Fig. 9, while the interaction vertices are given by Fig. 10.
z Y i s t H W b A j + 7 5 P / w s X O t m / 5 b L d x e F T F M c P W 2 D p r M p / t s U N 2 w k 5 Z i 3 F 2 y + 7 Z E 3 t 2 7 p w H 5 8 V 5 / W w d c 6 q Z V f a j n P c P q W C l 1 A = = < / l a t e x i t > . Feynman propagators, with the solid line being the φ cl φ cl 0 propagator, and the dashsolid line being the φ q φ cl 0 propagator.
t 9 C A J h B 4 h G d 4 h T f r y X q x 3 q 2 P 2 W j J K j L 7 8 A f W 5 w + q m 5 c F < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " q Z M S R n G X y d q X Q S L C e q j m v Y n 4 n 9 d N d X j l Z 0 w k q a a C z B a F K U c 6 R p M 6 U J 9 J S j Q f G 4 K J Z O a v i A y x a U S b 0 i q m B H f + 5 E X S O q u 5 h t + d V + v X R R 1 l O I Q j O A E X L q E O t 9 C A J h B 4 h G d 4 h T f r y X q x 3 q 2 P 2 W j J K j L 7 8 A f W 5 w + q m 5 c F < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " q Z M S R n G X y d q X Q S L C e q j m v Y n 4 n 9 d N d X j l Z 0 w k q a a C z B a F K U c 6 R p M 6 U J 9 J S j Q f G 4 K J Z O a v i A y x a U S b 0 i q m B H f + 5 E X S O q u 5 h t + d V + v X R R 1 l O I Q j O A E X L q E O t 9 C A J h B 4 h G d 4 h T f r y X q x 3 q 2 P 2 W j J K j L 7 8 A f W 5 w + q m 5 c F < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " q Z M S R n G X y d q X Q S L C e q j m v Y n 4 n 9 d N d X j l Z 0 w k q a a C z B a F K U c 6 R p M 6 U J 9 J S j Q f G 4 K J Z O a v i A y x a U S b 0 i q m B H f + 5 E X S O q u 5 h t + d V + v X R R 1 l O I Q j O A E X L q E O t 9 C A J h B 4 h G d 4 h T f r y X q x 3 q 2 P 2 W j J K j L 7 8 A f W 5 w + q m 5 c F < / l a t e x i t > i 24<l We now evaluate the loop correction to the advanced propagator, φ q φ cl , which we can see in terms of diagrams in Fig. 11.

3
Thermal initial density matrix for a free field 11 3.4 Initial density matrix for vacuum and n-particle states 14 3.5 Path integral with a free initial density matrix 15 3.6 Separating variables 17 3.7 One critical point for one initialization 17

Figure 3 .
Figure 3. Correlators for a single classical realisation (left) and averaged over initial conditions (right).

Figure 4 .
Figure 4.The variables to be integrated over on the real-time contour, after the initial conditions are fixed.

Figure 5 .
Figure 5.The classical correlator for a single initial condition, and the corresponding quantum averaged correlator.For λ = 0.0 (left) and 0.2 (right).

Figure 6 .
Figure 6.The full classical-statistical and quantum correlators (cl-cl) for a free and interacting theory at λ = 0.2.The figure on the right shows the result of subtracting the free propagator.The red line is the perturbative 1-loop result.

Figure 7 .
Figure 7. On the left, the full quantum correlators (cl-cl) for a free and interacting theory at λ = 4. On the right, when subtracting the free propagator.

Figure 8 .
Figure 8. Loop correction to the time-ordered two-point correlator, with the thick solid line being the Feynman propagator.
t e x i t s h a 1 _ b a s e 6 4 = " O G s b J c B j F 6 Y 5 I w I N l E U w F r 3 a O y 0 = " > A A A C M n i c b Z D L S s N A F I Y n 3 q 2 3 q k s 3 w S J U R E m K o M u i G 9 0 p W B W a E C b T k 3 Z w M h N m T o Q S 8 k x u f B L B h S 4 U c e t D O L 0 I 3 n 4 Y + P n O O Z w 5 f 5 w J b t D z n p y J y a n p m d m 5 + c r C 4 t L y y w a Z a B 5 t 4 f b I x I w Z e q B S q F L o 6 y N k b + L US P c L o s v V k b V m r f n D e X + N f 7 Y 1 M h Y Z 1 H 1 I e g o l q c g k Q l q T N v 3 M g w L q p E z A W U l y A 1 k l N 3 Q L r S t l T Q F E x b D k 0 t 3 y 5 K O m y h t n 0 R 3 S L 9 P F D Q 1 p p / G t j O l 2 D O / a w P 4 X 6 2 d Y 3 I Y F l x m O Y J k o 0 V J L l x U 7 i A / t 8 M 1 M B R9 a y j T 3 P 7 V Z T 1 q I 0 K b c s W G 4 P 8 + + a + 5 b O z 5 1 p / v 1 5 p H 4 z j m y A b Z J H X i k w P S J C f k j L Q I I 3 f k k b y Q V + f e e X b e n P d R 6 4 Q z n l k n P + R 8 f A L + o q s z < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " O G s b J c B j F 6 Y 5 I w I N l E U w F r 3 a O y 0 = " > A A A C M n i c b Z D L S s N A F I Y n 3 q 2 3 q k s 3 w S J U R E m K o M u i G 9 0 p W B W a E C b T k 3 Z w M h N m T o Q S 8 k x u f B L B h S 4 U c e t D O L 0 I 3 n 4 Y + P n O O Z w 5 f 5 w J b t D z n p y J y a n p m d m 5 + c r C 4 t L y y w a Z a B 5 t 4 f b I x I w Z e q B S q F L o 6 y N k b + L US P c L o s v V k b V m r f n D e X + N f 7 Y 1 M h Y Z 1 H 1 I e g o l q c g k Q l q T N v 3 M g w L q p E z A W U l y A 1 k l N 3 Q L r S t l T Q F E x b D k 0 t 3 y 5 K O m y h t n 0 R 3 S L 9 P F D Q 1 p p / G t j O l 2 D O / a w P 4 X 6 2 d Y 3 I Y F l x m O Y J k o 0 V J L l x U 7 i A / t 8 M 1 M B R9 a y j T 3 P 7 V Z T 1 q I 0 K b c s W G 4 P 8 + + a + 5 b O z 5 1 p / v 1 5 p H 4 z j m y A b Z J H X i k w P S J C f k j L Q I I 3 f k k b y Q V + f e e X b e n P d R 6 4 Q z n l k n P + R 8 f A L + o q s z < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " O G s b J c B j F 6 Y 5 I w I N l E U w F r 3 a O y 0 = " > A A A C M n i c b Z D L S s N A F I Y n 3 q 2 3 q k s 3 w S J U R E m K o M u i G 9 0 p W B W a E C b T k 3 Z w M h N m T o Q S 8 k x u f B L B h S 4 U c e t D O L 0 I 3 n 4 Y + P n O O Z w 5 f 5 w J b t D z n p y J y a n p m d m 5 + c r C 4 t L y y w a Z a B 5 t 4 f b I x I w Z e q B S q F L o 6 y N k b + L US P c L o s v V k b V m r f n D e X + N f 7 Y 1 M h Y Z 1 H 1 I e g o l q c g k Q l q T N v 3 M g w L q p E z A W U l y A 1 k l N 3 Q L r S t l T Q F E x b D k 0 t 3 y 5 K O m y h t n 0 R 3 S L 9 P F D Q 1 p p / G t j O l 2 D O / a w P 4 X 6 2 d Y 3 I Y F l x m O Y J k o 0 V J L l x U 7 i A / t 8 M 1 M B R9 a y j T 3 P 7 V Z T 1 q I 0 K b c s W G 4 P 8 + + a + 5 b O z 5 1 p / v 1 5 p H 4 z j m y A b Z J H X i k w P S J C f k j L Q I I 3 f k k b y Q V + f e e X b e n P d R 6 4 Q z n l k n P + R 8 f A L + o q s z < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " O G s b J c B j F 6 Y 5 I w I N l E U w F r 3 a O y 0 = " > A A A C M n i c b Z D L S s N A F I Y n 3 q 2 3 q k s 3 w S J U R E m K o M u i G 9 0 p W B W a E C b T k 3 Z w M h N m T o Q S 8 k x u f B L B h S 4 U c e t D O L 0 I 3 n 4 Y + P n O O Z w 5 f 5 w J b t D z n p y J y a n p m d m 5 + c r C 4 t L y y w a Z a B 5 t 4 f b I x I w Z e q B S q F L o 6 y N k b + L US P c L o s v V k b V m r f n D e X + N f 7 Y 1 M h Y Z 1 H 1 I e g o l q c g k Q l q T N v 3 M g w L q p E z A W U l y A 1 k l N 3 Q L r S t l T Q F E x b D k 0 t 3 y 5 K O m y h t n 0 R 3 S L 9 P F D Q 1 p p / G t j O l 2 D O / a w P 4 X 6 2 d Y 3 I Y F l x m O Y J k o 0 V J L l x U 7 i A / t 8 M 1 M B R9 a y j T 3 P 7 V Z T 1 q I 0 K b c s W G 4 P 8 + + a + 5 b O z 5 1 p / v 1 5 p H 4 z j m y A b Z J H X i k w P S J C f k j L Q I I 3 f k k b y Q V + f e e X b e n P d R 6 4 Q z n l k n P + R 8 f A L + o q s z < / l a t e x i t > i~✓(t 2 t 1 ) sin(! p [t 2 t 1 ]) !p < l a t e x i t s h a 1 _ b a s e 6 4 = " B l p h V t 2 E O B S v Y T j Q P c u S A q i b 8 o 8 = " > A A A C J 3 i c b Z B N S 8 N A E I Y 3 f l u / q h 6 9 B I t Q D 0 o i g p 5 E 9 O J R w a r Q h D D Z T t r F z S b s T o Q S + m + 8 + F e 8 C C q i R / + J 2 x r B r 4 G F h / e d Y X b e O J f C k O e 9 O W P j E 5 N T 0 z O z t b n 5 h c W l + v L K h c k K z b H F M 5 n p q x g M S q G w R Y I k X u U a I Y 0 l X s b X x 0 P / 8 g a 1 E Z k 6 p 3 6 O Y h O b p z 4 G 5 Y p e M m m b Z P k T t S v 0 + U k B r T T 2 P b m Q L 1 z G 9 v K P 7 n t Q t K 9 s N S q L w g V P x z U V J I l z J 3 G J r b E R o 5 y b 4 F 4 F r Y v 7 q 8 B z Y i s t H W b A j + 7 5 P / w s X O t m / 5 b L d x e F T F M c P W 2 D p r M p / t s U N 2 w k 5 Z i 3 F 2 y + 7 Z E 3 t 2 7 p w H 5 8 V 5 / W w d c 6 q Z V f a j n P c P q W i l 1 A = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " B l p h V t 2 E O B S v Y T j Q P c u S A q i b 8 o 8 = " > A A A C J 3 i c b Z B N S 8 N A E I Y 3 f l u / q h 6 9 B I t Q D 0 o i g p 5 E 9 O J R w a r Q h D D Z T t r F z S b s T o Q S + m + 8 + F e 8 C C q i R / + J 2 x r B r 4 G F h / e d Y X b e O J f C k O e 9 O W P j E 5 N T 0 z O z t b n 5 h c W l + v L K h c k K z b H F M 5 n p q x g M S q G w R Y I k X u U a I Y 0 l X s b X x 0 P / 8 g a 1 E Z k 6 p 3 6 O Y h O b p z 4 G 5 Y p e M m m b Z P k T t S v 0 + U k B r T T 2 P b m Q L 1 z G 9 v K P 7 n t Q t K 9 s N S q L w g V P x z U V J I l z J 3 G J r b E R o 5 y b 4 F 4 F r Y v 7 q 8 B z Y i s t H W b A j + 7 5 P / w s X O t m / 5 b L d x e F T F M c P W 2 D p r M p / t s U N 2 w k 5 Z i 3 F 2 y + 7 Z E 3 t 2 7 p w H 5 8 V 5 / W w d c 6 q Z V f a j n P c P q W i l 1 A = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " B l p h V t 2 E O B S v Y T j Q P c u S A q i b 8 o 8 = " > A A A C J 3 i c b Z B N S 8 N A E I Y 3 f l u / q h 6 9 B I t Q D 0 o i g p 5 E 9 O J R w a r Q h D D Z T t r F z S b s T o Q S + m + 8 + F e 8 C C q i R / + J 2 x r B r 4 G F h / e d Y X b e O J f C k O e 9 O W P j E 5 N T 0 z O z t b n 5 h h O b p z 4 G 5 Y p e M m m b Z P k T t S v 0 + U k B r T T 2 P b m Q L 1 z G 9 v K P 7 n t Q t K 9 s N S q L w g V P x z U V J I l z J 3 G J r b E R o 5 y b 4 F 4 F r Y v 7 q 8 B z Y i s t H W b A j + 7 5 P / w s X O t m / 5 b L d x e F T F M c P W 2 D p r M p / t s U N 2 w k 5 Z i 3 F 2 y + 7 Z E 3 t 2 7 p w H 5 8 V 5 / W w d c 6 q Z V f a j n P c P q W i l 1 A = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " B l p h V t 2 E O B S v Y T j Q P c u S A q i b 8 o 8 = " > A A A C J 3 i c b Z B N S 8 N A E I Y 3 f l u / q h 6 9 B I t Q D 0 o i g p 5 E 9 O J R w a r Q h D D Z T t r F z S b s T o Q S + m + 8 + F e 8 C C q i R / + J 2 x r B r 4 G F h / e d Y X b e O J f C k O e 9 O W P j E 5 N T 0 z O z t b n 5 h h O b p z 4 G 5 Y p e M m m b ZP k T t S v 0 + U k B r T T 2 P b m Q L 1 z G 9 v K P 7 n t Q t K 9 s N S q L w g V P x z U V J I l z J 3 G J r b E R o 5 y b 4 F 4 F r Y v 7 q 8 Bz Y i s t H W b A j + 7 5 P / w s X O t m / 5 b L d x e F T F M c P W 2 D p r M p / t s U N 2 w k 5 Z i 3 F 2 y + 7 Z E 3 t 2 7 p w H 5 8 V 5 / W w d c 6 q Z V f a j n P c P q W i l 1 A = = < / l a t e x i t > H W b A j + 7 5 P / w s X O t m / 5 b L d x e F T F M c P W 2 D p r M p / t s U N 2 w k 5 Z i 3 F 2 y + 7 Z E 3 t 2 7 p w H 5 8 V 5 / W w d c 6 q Z V f a j n P c P q W C l 1 A = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " C D Y j 3 r G j d 8 K z U u V 4 D n 3 0 e Y 0 p f / I = " > A A A C J 3 i c b Z B N S 8 N A E I Y 3 f l u / q h 6 9 B I t Q D 0 o i g p 5 E 9 O J R w a r Q h D D Z T t r F z S b s T o Q S + m + 8 + F e 8 C C q i R / + J 2 x r B r 4 G F h / e d Y X b e O J f C k O e 9 O W P j E 5 N T 0 z O z t b n 5 h H W b A j + 7 5 P / w s X O t m / 5 b L d x e F T F M c P W 2 D p r M p / t s U N 2 w k 5 Z i 3 F 2 y + 7 Z E 3 t 2 7 p w H 5 8 V 5 / W w d c 6 q Z V f a j n P c P q W C l 1 A = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " C D Y j 3 r G j d 8 K z U u V 4 D n 3 0 e Y 0 p f / I = " > A A A C J 3 i c b Z B N S 8 N A E I Y 3 f l u / q h 6 9 B I t Q D 0 o i g p 5 E 9 O J R w a r Q h D D Z T t r F z S b s T o Q S + m + 8 + F e 8 C C q i R / + J 2 x r B r 4 G F h / e d Y X b e O J f C k O e 9 O W P j E 5 N T 0 z O z t b n 5 h H W b A j + 7 5 P / w s X O t m / 5 b L d x e F T F M c P W 2 D p r M p / t s U N 2 w k 5 Z i 3 F 2 y + 7 Z E 3 t 2 7 p w H 5 8 V 5 / W w d c 6 q Z V f a j n P c P q W C l 1 A = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " C D Y j 3 r G j d 8 K z U u V 4 D n 3 0 e Y 0 p f / I = " > A A A C J 3 i c b Z B N S 8 N A E I Y 3 f l u / q h 6 9 B I t Q D 0 o i g p 5 E 9 O J R w a r Q h D D Z T t r F z S b s T o Q S + m + 8 + F e 8 C C q i R / + J 2 x r B r 4 G F h / e d Y X b e O J f C k O e 9 O W P j E 5 N T 0 z O z t b n 5 h 2 M A v l C J 8 K X d S I = " > A A A C A X i c b V D L S s N A F L 2 p r 1 p f U T e C m 8 E i u L E k I u q y 6 M Z l B f u A J p T J d N I O n U z C z E Q o I W 7 8 F T c u F H H r X 7 j z b 5 y 2 W W j r g Y H D O f d w 5 5 4 g 4 U x p x / m 2 S k v L K 6 t r 5 f X K x u b W 9 o 6 9 u 9 d S c S o J b Z K Y x 7 I T Y E U 5 E 7 S p m e a 0 k 0 i K o 4 D T d j C 6 m f j t B y o V i 8 W 9 H i f U j / B A s J A R r I 3 U s w 9 O v V B i k j G P m 1 A f 5 9 m F N w y w z H t 2 1 a k 5 U 6 B F 4 h a k C g U a P f v L 6 8 c k j a j Q h G O l u q 6 T a D / D U j P C a V 7 x U k U T T E Z 4 Q L u G C h x R 5 W f T C 3 J 0 b J Q + C m N p n t B o q v 5 O Z D h S a h w F Z j 2 M A v l C J 8 K X d S I = " > A A A C A X i c b V D L S s N A F L 2 p r 1 p f U T e C m 8 E i u L E k I u q y 6 M Z l B f u A J p T J d N I O n U z C z E Q o I W 7 8 F T c u F H H r X 7 j z b 5 y 2 W W j r g Y H D O f d w 5 5 4 g 4 U x p x / m 2 S k v L K 6 t r 5 f X K x u b W 9 o 6 9 u 9 d S c S o J b Z K Y x 7 I T Y E U 5 E 7 S p m e a 0 k 0 i K o 4 D T d j C 6 m f j t B y o V i 8 W 9 H i f U j / B A s J A R r I 3 U s w 9 O v V B i k j G P m 1 A f 5 9 m F N w y w z H t 2 1 a k 5 U 6 B F 4 h a k C g U a P f v L 6 8 c k j a j Q h G O l u q 6 T a D / D U j P C a V 7 x U k U T T E Z 4 Q L u G C h x R 5 W f T C 3 J 0 b J Q + C m N p n t B o q v 5 O Z D h S a h w F Z j 2 M A v l C J 8 K X d S I = " > A A A C A X i c b V D L S s N A F L 2 p r 1 p f U T e C m 8 E i u L E k I u q y 6 M Z l B f u A J p T J d N I O n U z C z E Q o I W 7 8 F T c u F H H r X 7 j z b 5 y 2 W W j r g Y H D O f d w 5 5 4 g 4 U x p x / m 2 S k v L K 6 t r 5 f X K x u b W 9 o 6 9 u 9 d S c S o J b Z K Y x 7 I T Y E U 5 E 7 S p m e a 0 k 0 i K o 4 D T d j C 6 m f j t B y o V i 8 W 9 H i f U j / B A s J A R r I 3 U s w 9 O v V B i k j G P m 1 A f 5 9 m F N w y w z H t 2 1 a k 5 U 6 B F 4 h a k C g U a P f v L 6 8 c k j a j Q h G O l u q 6 T a D / D U j P C a V 7 x U k U T T E Z 4 Q L u G C h x R 5 W f T C 3 J 0 b J Q + C m N p n t B o q v 5 O Z D h S a h w F Z j

Figure 10 .
Figure 10.Feynman diagrams for the interactions, with the solid line representing φ cl , and the dashed line corresponding φ q .