Tree-level correlations in the strong field regime

We consider the correlation function of an arbitrary number of local observables in quantum field theory, in situations where the field amplitude is large. Using a quasi-classical approximation (valid for a highly occupied initial mixed state, or for a coherent initial state if the classical dynamics has instabilities), we show that at tree level these correlations are dominated by fluctuations at the initial time. We obtain a general expression of the correlation functions in terms of the classical solution of the field equation of motion and its derivatives with respect to its initial conditions, that can be arranged graphically as the sum of labeled trees where the nodes are the individual observables, and the links are pairs of derivatives acting on them. For 3-point (and higher) correlation functions, there are additional tree-level terms beyond the quasi-classical approximation, generated by fluctuations in the bulk.


JHEP09(2017)055 1 Introduction
A common question in many areas is the evaluation of correlations between several measurements, given the microscopic dynamics of the system. This problem comes up for instance in the calculation of cosmological perturbations [1,2], or in nuclear physics, e.g. in heavy ion collisions [3,4]. These situations have in common that they are described by some underlying quantum field theory, and that the system starts from some supposedly known initial state. Then, it evolves under the effect of the self-interactions of the fields, and possibly the couplings to some external sources. Thereafter, we define an observable as some local operator O(φ(x)) constructed from the fields of the theory. Our goal in this paper is to study (at leading order in the couplings) the correlations between measurements of this observable at several space-time points x 1 , x 2 , · · · , x n , (1.1) In this definition, we have assumed that the system is initially in the perturbative vacuum state (see the appendices B and C for a discussion of the changes with other types of initial states). The subscript "c" indicates the connected part of this expectation value, i.e. the terms that contain the actual correlations. An important restriction in our discussion will be to consider only points x i that all have space-like separations, i.e. (x i − x j ) 2 < 0. Physically, this means that the measurement performed at one of the points cannot influence the outcome of the measurement at another point. Therefore, the correlations between the measurements are entirely due to the past evolution of the system. The main difference compared to the more familiar problem of computing cross-sections is that the final state is not prescribed. Instead, one wishes to calculate the expectation value of some operators, by summing over all the possible final states given a prescribed initial state. 1 One could in principle perform such a calculation by using the usual techniques for calculating transition amplitudes, by first writing (1.2) (The double sum over final states could be reduced to a single sum if we use eigenstates of the observable O.) All the steps involved in this pedestrian approach can be encapsulated in a set of diagrammatic rules known as the Schwinger-Keldysh formalism, or in-in formalism [5][6][7][8] (see also [9,10]), that roughly consists in two copies of the usual Feynman rules (one that gives the factor f in , and one -complex conjugated-that gives the factor in f ), plus some additional rules that give the final state sum. The approach we follow in this paper bears heavily on earlier works in the field of heavy ion collisions [11][12][13][14] and in cosmology [2,15].
Although the approach used in this paper and the final result are quite general, the intermediate derivations are a bit cumbersome. In order to keep the notations as light as JHEP09(2017)055 possible, we consider a theory with a single real scalar field φ, whose Lagrangian density is given by , (1.3) where V (φ) contains the self-interactions of the field and J(x) is an external source. The quantum field theory described by the Lagrangian density (1.3) has a well known perturbative expansion. However, we are interested in this paper in the strong field regime, where this perturbative expansion is insufficient. This corresponds to fields for which the interaction term V (φ) is as large as the kinetic term (∂ µ φ)(∂ µ φ), implying that the interactions cannot be treated as a perturbation. For instance, for a theory with a quartic interaction term V (φ) ∼ g 2 φ 4 , this occurs for fields φ ∼ Q g , (1.5) where Q is some typical momentum scale in the problem (the field φ has the dimension of a momentum in four space-time dimensions). For fields driven by an external source, the term Jφ must also have the same order of magnitude as the kinetic and self-interaction terms, and thus the order of magnitude of the source that may create this large field 2 is (1. 6) In order to see where the usual perturbative expansion fails, let us recall the power counting in the case of a theory with a g 2 φ 4 interaction term. The order of magnitude of a connected graph G with n E external legs, n L loops and n J external sources is given by (1.7) Eq. (1.6) implies that gJ cannot be treated as a small expansion parameter. Any quantity must therefore be evaluated non-perturbatively in the number of sources n J . In contrast, the hierarchy of the contributions based on the number of loops in the Feynman diagrams survives. Our aim is to calculate the correlation function (1.1) to lowest order in the number of loops (i.e. at tree level) but to all order in n J (or equivalently to all orders in the field amplitude). With large fields, the correlation function (1.1) is suppressed compared to the uncorrelated part, because each propagator connecting a pair of observables costs a factor g 2 (the endpoints of such a link replace two fields of order g −1 ). In order to fully connect the n observables, n − 1 propagators are necessary, suppressing the correlation by a factor g 2(n−1) with respect to the uncorrelated part. Thus, in a certain sense, the correlation JHEP09(2017)055 function may be viewed as a higher order correction. At leading order, 1-point functions (i.e. the expectation value of a local observable) are obtained from a classical solution of the field equation of motion, obeying a retarded boundary condition that depends on the initial state (φ =φ = 0 when the initial state is the perturbative vacuum) [11]. Their nextto-leading order (NLO) corrections are also known [13], and can be expressed in terms of functional derivatives of this classical field with respect to its initial condition. A similar result was proven in [14] for the NLO correction to a 2-point function, that contains its tree-level correlated part.
In this paper, we assess to what extent a similar representation exists for the correlated part of higher n-point functions at tree-level. We prove that, within the quasi-classical approximation (defined in the section 5), the tree-level correlation between any number of arbitrary local observables takes the following graphical form: where the nodes are the observables evaluated at a classical solution of the field equation of motion, and the links are pairs of derivatives with respect to the initial condition of the classical field. The paper is organized as follows. In the section 2, we couple the observable O(φ(x)) to a fictitious source z(x) in order to define a generating functional whose derivatives give the expectation for the combined measurement of O(φ(x)) at several points. We explain how to obtain this generating functional diagrammatically by extending the Schwinger-Keldysh formalism with an extra vertex that corresponds to the observable, and then we rephrase this diagrammatic expansion in the retarded-advanced basis, that will be more convenient later. In the section 3, we study the first derivative of the generating functional. At leading order, it can be expressed in terms of a pair of fields that obey classical equations of motion and satisfy non-trivial z-dependent boundary conditions. In the following section 4, we set up an expansion of these fields order by order in powers of the function z(x), and we determine explicitly the first two orders, that give the 1-point and 2-point correlation functions. We also show how to rewrite the 2-point function in terms of derivatives of a classical field with respect to its initial condition. In the section 5, we consider the quasi-classical approximation, and we determine the solution to all orders in z in this approximation. We show that the tree level n-point correlation function in this approximation can be represented as the sum of all trees with n labeled nodes (the n instances of the observable), and pairwise links that are differential operators with respect to the initial condition of the classical field. The section 6 illustrates on the example of the 3-point function the type of contributions that arise at tree level beyond the quasi-classical approximation. Summary and conclusions are in the section 7. In the appendix A, we derive some technical points used in the main part of the paper, while in the appendices B and C we discuss other types of initial states.

Definition
For simplicity, we will take the points x i where the measurements are performed to lie on the same surface of constant time x 0 = t f , but the final results would be valid for any locally space-like surface (in order to ensure that there is no causal relation between the points x i ).
One can encapsulate all the correlation functions (1.1) into a generating functional defined as follows: 3 where the argument of the field in O is x ≡ (t f , x). From this generating functional, the correlation functions are obtained by differentiating with respect to z(x 1 ), · · · , z(x n ) and by setting z ≡ 0 afterwards. In order to remove the uncorrelated part of the n-point function, we should differentiate the logarithm of F , i.e.
Note that since the final surface is space-like and the operators O(φ(x)) are local, there is no need for a specific ordering of the exponential in eq. (2.1). The observable O(φ(x)) is made of the field in the Heisenberg picture, φ(x), that can be related to the field φ in (x) of the interaction picture as follows: where U (t 1 , t 2 ) is an evolution operator given in terms of the interactions by the following formula We can therefore rewrite the generating functional solely in terms of the interaction picture field φ in , where P denotes a path ordering on the following time contour: Figure 1. Diagrammatic rules for the extended Schwinger-Keldysh formalism that gives the generating functional. The Feynman rules shown here for the self-interactions correspond to a g 2 φ 4 /4! interaction term. In this illustration, we have assumed that the observable is quartic in the field when drawing the corresponding vertex (proportional to z(x)).
P orders the operators from left to right starting from the end of the contour. We denote by φ in+ the field that lives on the upper branch and by φ in− the field on the lower branch (the minus sign in front of the term L int (φ in− (x)) comes from the fact that the lower branch is oriented from +∞ to −∞). The operator O(φ in (x)) lives at the final time of this contour, and could either be viewed as made of fields of type + or of type − (the two choices lead to the same results). Note that in eq. (2.5), the spacetime integration is a priori extended to the domain located below the final surface (the observable O therefore lives on the upper time boundary of the integration domain). However, this restriction is not really necessary: by causality, the contribution of the domain located above t f would cancel anyhow.

Expression in the Schwinger-Keldysh formalism
Since the initial state is the vacuum, the generating functional defined in eq. (2.5) can be represented diagrammatically as the sum of all the vacuum-to-vacuum graphs (i.e. graphs without external legs) in the Schwinger-Keldysh formalism (also known as the in-in formalism), extended by an extra vertex that corresponds to the insertions of the observable O. Let us recall here that the Schwinger-Keldysh diagrammatic rules consist in having two types of interaction vertices (+ and − depending on which branch of the contour the vertex lies on, the − vertex being the opposite of the + one) and four types of bare propagators (G 0 ++ , G 0 −− , G 0 +− and G 0 −+ ) depending on the location of the endpoints on the contour. The additional vertex exists only on the final surface, at the time t f . It is accompanied by a factor z(x), and has as many legs as there are fields in O. There is only one kind of this vertex (we can decide to call it + or − without affecting anything). We recapitulate these Feynman rules in the figure 1. In the case of the vacuum initial state, the propagators have the following explicit expressions:

JHEP09(2017)055
where we have used the following compact notation Note that when we set z ≡ 0, these diagrammatic rules fall back to the pure Schwinger-Keldysh formalism, for which it is known that all the connected vacuum-to-vacuum graphs are zero. This implies that in accordance with the fact that this should be 0 in 0 in = 1.

Retarded-advanced representation
In order to clarify what approximations may be done in the large field regime, it is useful to use a different basis of fields by introducing [16,17] The half-sum φ 2 in a sense captures the classical content (plus some quantum corrections), while the difference φ 1 is purely quantum (because it represents the different histories of the fields in the amplitude and in the complex conjugated amplitude). To see how the Feynman rules are modified in terms of these new fields, let us start from where the matrix Ω reads: In terms of this matrix, the new propagators are obtained as follows Explicitly, they read Note that G 0 21 is the bare retarded propagator, while G 0 12 is the bare advanced propagator. The vertices in the new formalism (here written for a quartic interaction g 2 φ 4 /4!) are given by 14)

JHEP09(2017)055
where More explicitly, we have: The vertices not listed explicitly here are obtained by permutations.) Moreover, the rules for an external source in the retarded-advanced basis are: Finally, note that the observable depends only on the field φ 2 , i.e. O = O(φ 2 ). Indeed, the fields φ + and φ − represent the field in the amplitude and in the conjugated amplitude. Their difference φ 1 should vanish when a measurement is performed.

First derivative of ln F
Differentiating the generating functional with respect to z(x) amounts to exhibiting a vertex O at the point x at the final time (as opposed to weighting this vertex by z(x) and integrating over x). Furthermore, by considering the logarithm of the generating functional rather than F itself, we have only diagrams that are connected to the point x, as shown in this representation: where the gray blob is a sum of graphs constructed with the Feynman rules of the figure 1, or their analogue in the retarded-advanced formulation. Therefore, these graphs still depend implicitly on z. Note that this blob does not have to be connected.

Tree level expression
Without further specifying the content of the blob, eq. (3.1) is valid to all orders, both in z and in g. At lowest order in g (tree level), a considerable simplification happens because the blob must be a product of disconnected subgraphs, one for each line attached to the vertex O(φ(x)): where now each of the light colored blob is a connected tree 1-point diagram. In the retarded-advanced formalism, there are two of these 1-point functions, that we will denote JHEP09(2017)055 φ 1 and φ 2 . At tree level, they can be defined recursively by the following pair of coupled integral equations: In these equations, O is the derivative of the observable with respect to the field, Ω is the space-time domain comprised between the initial and final times, and we denote For an interaction Lagrangian − g 2 4! φ 4 + Jφ, this difference reads In terms of these fields, we have i.e. simply the observable O evaluated on the field φ 2 (x) (but this field depends on z to all orders, via the boundary terms in eqs. (3.3)).

Classical equations of motion
Using the fact that G 0 12 and G 0 21 are Green's functions of + m 2 , obeying the following identities while G 0 22 vanishes when acted upon by this operator, we see that φ 1 and φ 2 obey the following classical field equations of motion: Note that here the point x is located in the "bulk" Ω; this is why the observable does not enter in these equations of motion. In fact, the observable enters only in the boundary conditions satisfied by these fields on the hypersurface x 0 = t f . For later reference, let us

JHEP09(2017)055
also rewrite these equations of motion in the specific case of a scalar field theory with a g 2 φ 4 /4! interaction term and an external source J:

Boundary conditions
The equations of motion (3.9) are easier to handle than the integral equations (3.3), but they must be supplemented with boundary conditions in order to define uniquely the solutions. The standard procedure for deriving the boundary conditions is to consider the combination G 0 12 (x, y) ( y + m 2 ) φ 1 (y), and let the operator y + m 2 act alternatively on the right and on the left, By subtracting these equations and integrating over y ∈ Ω, we obtain The second term of the right hand side is a total derivative thanks to Therefore, this term can be rewritten as a surface integral extended to the boundary of the domain Ω. With reasonable assumptions on the spatial localization of the source J(x) that drives the field, we may disregard the contribution from the boundary at spatial infinity. The remaining boundaries are at the initial time t i and final time t f , (3.14) Note that the boundary term vanishes at the initial time t i , because G 0 12 is the advanced propagator. Likewise, we obtain the following equation for φ 2 : Then, the boundary conditions at t i and t f are obtained by comparing eqs. (3.3) and (3.14)-(3.15). At the final time t f , the boundary condition is

JHEP09(2017)055
At the initial time t i , we must have Some simple manipulations lead to the following equivalent form From the explicit form of the propagators G 0 +− and G 0 −+ (given in eqs. (2.6)), we see that, at the initial time, φ 2 + 1 2 φ 1 has no positive frequency components, and φ 2 − 1 2 φ 1 has no negative frequency components. An equivalent way to state this boundary condition is in terms of the Fourier coefficients of the fields φ 1,2 . Let us decompose them at the time t i as follows, In terms of the coefficients introduced in this decomposition, the boundary conditions at the initial time read: The equations of motion (3.9), accompanied by the boundary conditions (3.16) and (3.20), are equivalent to the formulation of [15] (this reference uses the fields φ ± of the in-in formalism, instead of the fields φ 1,2 of the retarded-advanced representation that we are using here).

Setup of the expansion
In the tree level approximation, the fields φ 1,2 obey classical equations of motion and satisfy coupled boundary conditions (both at the initial and final times), a problem which is usually difficult to solve, even numerically. Moreover, one should keep in mind that the z-dependence of the solutions arises entirely from the boundary condition at the final time.
A first approach, that we shall pursue in this section, it to expand the solutions in powers of z, by writing them as follows: 2 (x; x 1 , x 2 ) + · · · (4.1)

JHEP09(2017)055
By construction, the coefficients of this expansion are symmetric functions of the x i 's, and they are nothing but the functional derivatives of φ 1,2 (x) with respect to z, The n-point correlation function is obtained by starting from the first derivative of ln F , equal to O(φ 2 ), and by differentiating it n−1 times with respect to z. Using Faà di Bruno's formula, this leads to where the notations used in this equation are the following: • π denotes one of these partitions, and |π| is its cardinal, i.e. the number of blocks into which {2 · · · n} is partitioned, • σ denotes a block in the partition π, and |σ| the number of elements in this block.
From eq. (4.3), it is clear that all the correlation functions can be obtained from the coefficients in the expansion of φ 2 in powers of z.
The coefficients φ (n) 1,2 can be determined iteratively as follows. By inserting eqs. (4.1) in the equations of motion (3.9) (or eqs. (3.10) in the specific case of a φ 4 interaction), and by using the fact that the functions 1, z(x 1 ), z(x 1 )z(x 2 ), · · · are linearly independent, we obtain a system of coupled equations for the coefficients. In a similar fashion, we obtain boundary conditions at the initial and final time for the coefficients. These equations form a "triangular" hierarchy, starting with equations that involve only φ 1,2 with n < n. In the rest of this section, we solve the first two orders of this hierarchy, in order to obtain the 1-point and 2-point correlation functions.

Order zero
The order 0 in z is obtained by setting z ≡ 0 in the equations of motion (3.9) and in the boundary conditions (3.16)- (3.18). This leads to considerable simplifications. Firstly, the boundary condition at t f for the coefficient φ Then, the boundary conditions at t i tell us that i.e. the usual Euler-Lagrange equation for the theory under consideration. Therefore, φ (0) 1 is zero, and φ 2 is the solution of the classical equation of motion that vanishes (as well as its first time derivative) at the initial time. The expectation value of the observable at tree level is simply given by (4.8) Since this object will appear repeatedly in the following, let us introduce a compact graphical representation: The shaded blob in fact encapsulates an infinite series of tree diagrams, a glimpse of which is given by the second equality.

Order one 4.3.1 Equations of motion and boundary conditions
Let us now consider the next order in z. Firstly, we obtain the following two equations of motion, where we have used the fact that φ while those at the initial time are JHEP09(2017)055

Solution in terms of mode functions
The solution of this problem can be constructed as follows. Firstly, let us introduce a complete basis of solutions of the linear equation of motion that appears in eq. (4.10). A convenient choice will be to choose the functions of this basis such that they coincide with plane waves at the initial time, These functions are sometimes called mode functions. The boundary condition at x 0 → t i means that both the value of the function and that of its first time derivative coincide with those of the indicated plane wave. Thus, the functions a +k (x) contain only positive frequency modes at t i , and the functions a −k have only negative frequencies. 4 The solution for φ 1 can be expressed as a linear combination of these mode functions as follows,  Since the mode functions a +k (x) and a −k (x) have been defined in such a way that they contain only positive or negative frequencies at t i , respectively, it is immediate to write the solution for φ Note that the underlined combination, that we denote G 22 , is nothing but the symmetric 22 propagator dressed by the background field φ 2 . The 2-point correlation function is then given by (4.16) Note that, although the final expression for the correlation function is symmetric (as expected since the two observables with a space-like separation commute), the intermediate calculations break this manifest symmetry. Eq. (4.16) generalizes to the non-linear strong field regime the "equation of motion method" (see [18] for instance) that has been devised for the calculation of 2-point correlations of cosmological density fluctuations.

Expression in terms of derivatives w.r.t. the initial field
Eq. (4.16) can also be rewritten in a different way, that we will generalize to the case of the n-point correlation in the quasi-classical approximation. Firstly, let us generalize the classical field φ (0) 2 (x) (that has null initial conditions) into a classical field Φ obeying the same equation of motion but with generic initial condition Φ ini at t i . Such a field obeys the following Green's formula, where G 0 R ≡ G 0 ++ − G 0 +− is the bare retarded propagator. Note that, by definition Then, let us write a similar integral representation for the mode functions a ±k (x) introduced above, By comparing these two integral representations, we see that we can write the following formal relationship between Φ and a ±k : . (4.20) In other words, the mode function a ±k (x) can be viewed as the functional derivative of the retarded classical field Φ with respect to its initial condition, reweighted by the plane wave exp(∓ik · x). This relationship implies that (4.21) Thus, the tree level 2-point correlation function can be written as . (4.22)

JHEP09(2017)055
The ⊗ product 5 introduced here is just a compact notation for the symmetrized integration over the momentum k. In this formulation, the 2-point correlation function is obtained by linking two copies of the observable (at the points x 1 and x 2 ), with a "link operator" made of a pair of derivatives with respect to the initial condition of the classical field Φ, one acting on each end of the link. We will represent diagrammatically this operation as follows: It is important to keep in mind that this compact representation is not a Feynman diagram, but rather a shorthand for an infinite series of tree Feynman diagrams, made of two copies of the graphs that appear in eq. (4.9) connected by a bare G 0 22 propagator in all the possible ways: This representation of the tree-level 2-point correlation function was already contained in the NLO result obtained in [14], but the retarded-advanced basis used in the present paper simplifies considerably its derivation.
The main result of this section is that all these graphs can be obtained from the 1point function evaluated for a classical field Φ with generic initial condition Φ ini , by taking functional derivatives with respect to this initial condition. The "link" that connects the two handles created by these derivatives encodes the zero-point fluctuations in the initial vacuum state. For higher-point correlation functions at tree level, this result is not true in general. Starting with the 3-point correlation function (see the section 6), there are additional contributions that cannot be expressed in terms of functional derivatives with respect to the initial classical field Φ ini . However, as we shall see in the next section, these additional contributions are absent in the quasi-classical approximation.

Quasi-classical approximation
Until now, our counting was simply based on the fact that a large external source J leads to large fields φ ± , but no additional approximation was made in the calculation of the 1-point and 2-point correlation functions at tree level in the previous section. Instead of pursuing the very cumbersome expansion in powers of z that we have used so far, we 5 From its definition, the ⊗ product is symmetric,

JHEP09(2017)055
consider in this section an approximation that allows a compact formal solution to all orders in z. In this subsection, we give only a very sketchy motivation for this approximation, and then go straight to its implications for the calculation of correlations. Afterwards, a lengthier discussion of its validity will be provided later in this section (after we have derived expressions for the fields φ 1 and φ 2 , that are necessary for this). Let us first recall that the fields φ + and φ − represent, respectively, the space-time evolution of the field in amplitudes and in conjugate amplitudes. The fact that they are distinct leads to quantum interferences when squaring amplitudes, an effect controlled by . When fields are strong (e.g., because of a strong external source, or because the initial condition has a large coherent field), we may intuitively expect the system to behave quasiclassically, i.e. that the quantum interferences are a small correction. In more technical terms, this translates into a small φ 1 ≡ φ + − φ − compared to φ ± . But it turns out that the situation is more complicated, because φ 1,2 also have an implicit dependence on the dummy source z(x). For this to lead to useful simplifications, we must request that In this regime, that we will call the quasi-classical approximation (QCA), we can approximate the equations of motion 6 (3.9) by keeping only the lowest order in φ 1 . This amounts to keeping only the terms linear in φ 1 in eq. (3.4) (in the case of a φ 4 theory, it means dropping the φ 3 1 φ 2 term in eq. (3.5)). In the approximation, they read while the boundary conditions are still given by (3.16) and (3.18). The problem one must now solve is illustrated in the figure 2. The field φ 1 obeys a linear equation of motion (dressed by the field φ 2 , although this aspect is not visible in the figure), with an advanced boundary condition that depends on φ 2 . In parallel, the field φ 2 obeys a non-linear equation of motion, with a retarded boundary condition that depends on φ 1 . As we shall show in the next subsection, this tightly constrained problem admits a formal solution, valid to all orders in the function z, in the form of an implicit functional equation for the first derivative of ln F [z].

Formal solution
The equation of motion for φ 1 (first of eqs. (5.3)) is formally identical to the first of eqs. (4.10). We can therefore mimic eq. (4.14) and write directly φ 1 as follows: Figure 2. Relationship between the fields φ 1 and φ 2 in the quasi-classical approximation.
where the mode functions a ±k should now be defined with φ 2 as the background, rather than φ 2 . To obtain eq. (5.4), it was crucial to have a linear equation of motion for φ 1 , a consequence of the quasi-classical approximation. The above equation formally defines φ 1 (x) in the bulk, x ∈ Ω, in terms of the field φ 2 at the final time. It is important to note that this equation is valid to all orders in z, contrary to the equations encountered in the expansion in powers of z that we have used in the previous section. Besides the explicit factor z(u), the right hand side contains also an implicit z dependence (to all orders in z) in the field φ 2 (t f , u) and in the mode functions a ±k (since they evolve on top of the background φ 2 ).
Then, using the boundary condition at the initial time, we obtain the following expression for the field φ 2 at t i , where the arrows indicate on which side the T ±k operators act. This expression for the field φ 2 at the time t i can be used as initial condition for the first of eqs. (5.3). The next step is to note that the field φ 2 (x) that satisfies this equation of motion, and has the initial condition φ 2 (t i , y) is formally given by This formula follows from the fact that the derivative δ/δΦ ini is the generator for shifts of the initial condition of Φ; its exponential is therefore the corresponding translation operator. The same formula applies also to any function of the field, since the exponential operator shifts Φ ini in any instance of Φ on its right. In particular, we have . Substituting φ 2 (t i , y) by eq. (5.5) inside the exponential, this leads to the first derivative of ln F , we see that it obeys the following recursive formula Let us make a technical remark about the scope of the ← T±k derivative (contained in the ⊗ product) that acts on the factor D[u; z] inside the exponential. This derivative originates from the terms in a ±k (t f , u) O (φ 2 (t f , u)) in eq. (5.5). From this origin, it is clear that ← T±k must establish a link between a point on the initial surface and the point u on the final surface, as illustrated in the figure 3, where the black line is the function a ±k (t f , u)). As we shall see in the next section, the recursive expansion of eq. (5.10) leads to a tree structure where the nodes are observables O (Φ(t f , x)). The above discussion tells us that when performing this expansion, the ⊗ product should not be "distributed" to all the nodes inside D[u; z], but applied only to the node of coordinate u.

Realization of the quasi-classical approximation
Let us now return on the condition φ 1 φ 2 , that was used in the derivation of eq. (5.10), in order to see a posteriori when it is satisfied. Firstly, let us make a number of simple remarks:

JHEP09(2017)055
• The boundary condition at the initial time t i states that the Fourier coefficients of φ 1 and φ 2 are proportional. Therefore, at least in the immediate vicinity of the initial time, φ 1 and φ 2 must have comparable orders of magnitude, to all orders in z(x), • The large external source J enters only in the equation of motion for φ (0) 2 . It is neither present in the equation for φ 1 , nor in the equations that drive the higher order coefficients of φ 2 .
From these two facts, we conclude that after some amount of time evolution, φ (0) 2 becomes much larger than φ (0) 1 (this conclusion is correct), while we naively expect that φ for n ≥ 1 (as we shall see, instabilities change this conclusion).
For a more elaborate discussion, we can use eq. (5.4) for φ 1 . For φ 2 , the initial condition at t i is given by eq. (5.5). For the sake of this discussion, it is sufficient to use a sourceless linearized solution for φ 2 in the bulk, that reads u)) .
(5.11) (This expression captures correctly the order of magnitude of all φ (n) 2 , except for n = 0 since the source has been removed -but for n = 0, we already know that the quasiclassical approximation is valid.) First of all, a comparison between eqs. (5.4) and (5.11) indicates that φ 1 and φ 2 have the same order in the coupling constant g, since they are made of the same building blocks (the only difference is the sign between the two terms of the integrand, and an irrelevant overall factor 1 2 ). However, a hierarchy between φ 1 and φ 2 arises dynamically at every order in z(x) when the classical solutions of the field equation of motion (4.7) are unstable. Such instabilities are fairly generic in several quantum field theories; in particular the scalar field theory with a φ 4 coupling that we are using as example throughout this paper is known to have a parametric resonance [20,21]. Since the mode functions a ±k are linearized perturbations on top of the classical field φ 2 , an instability of the classical solution φ 2 is equivalent to the fact that some of the mode functions grow exponentially with time, as exp(µ(x 0 − t i )) (where µ is the Lyapunov exponent). Thus, since eq. (5.11) is bilinear in the mode functions, we expect that for sufficiently large t f − t i . In order to estimate the ratio φ 2 /φ 1 at intermediate times, one may use the following reasoning. The antisymmetric combination of mode functions that enters in eq. (5.4) is the advanced propagator G A in the background φ 2 . This advanced propagator may also JHEP09(2017)055 be expressed in terms of a different set of mode functions b ±k defined to be plane waves at the final time t f , (5.14) In terms of these alternate mode functions, we also have In the presence of instabilities, these backward evolving mode functions grow when x 0 decreases away from t f , as exp(µ(t f − x 0 )) (in this sketchy argument, the Lyapunov exponent µ is assumed here to be the same for the forward and backward mode functions). This implies 16) and the following magnitude for the ratio φ 2 /φ 1 at intermediate times Thus, with instabilities and non-zero Lyapunov exponents, the strong field approximation is generically satisfied thanks to the exponential growth of perturbations over the background.
In the appendix C, we will discuss another situation where this approximation is also satisfied, namely when the initial state is a mixed state with large occupation number.

Expansion of eq. (5.10) in powers of z
Although eq. (5.10) cannot be solved explicitly, it is fairly easy to obtain a diagrammatic representation of its solution. For this, let us introduce the following graphical notations: At the order 0 in z, we just need to set z ≡ 0 inside the exponential, to obtain Then, we proceed recursively. We insert the 0-th order result in the exponential, and expand to order 1 in z, leading to the following result at order 1: The next two iterations give: Note that the sum of the weights for the trees with n + 1 nodes (one of them being the root node x 1 ), This functional identity may be viewed as a structureless version of eq. (5.10), in which the function z(x) is replaced by a scalar variable z. Eq. (5.25) leads to (see [22], pp. 127-128) which is the number of trees with n + 1 labeled nodes (Cayley's formula) divided by the number of ways of reshuffling the n nodes that are integrated out.

Correlation functions
The n-point correlation function is obtained by differentiating this expression n − 1 times, with respect to z(x 2 ), · · · , z(x n ), and by setting z ≡ 0 afterwards. This selects all the JHEP09(2017)055 trees with n distinct labeled 7 nodes (including the node at x 1 ). Moreover, since derivatives commute, these successive differentiations eliminate the symmetry factors, leading to (5.27) The number of trees contributing to this sum is equal to n n−2 . Eq. (5.27) tells us that, at tree level in the quasi-classical regime, all the n-point correlation functions are entirely determined by the functional dependence of the solution of the classical field equation of motion with respect to its initial condition. Moreover, this formula provides a way to construct explicitly the correlation functions in terms of functional derivatives with respect to the initial field. In this tree representation, the number of links reaching a node is the number derivatives with respect to Φ ini that act on the corresponding O(Φ). When the observable is a composite operator in terms of the fields of the theory, this includes a large number of contributions. For instance, when O(Φ(x)) is cubic in the field Φ(x), a node with three links would contain the following terms: where the three dots inside the blob represent the three fields Φ it contains. Instead of the functional derivation of the representation (5.27) that we have performed in this section, one could in principle have iterated on the z-expansion introduced in the previous section (after simplifying the equations of motion based on φ 1 φ 2 ). This approach produces a sum of terms corresponding to the right hand side of eq. (5.28) that, after some hefty combinatorics, one may recognize as being the compact form of the left hand side of eq. (5.28).
In the quasi-classical approximation, the final state correlations are entirely due to quantum fluctuations in the initial state, that are encoded in the function G 0 22 (x, y). If the initial state is the vacuum, as we have assumed in this paper, it reads The support of this function is dominated by distances |x − y| smaller than the Compton wavelength m −1 . Thus, in the tree representation of eq. (5.27), a link between the points x i and x j is nonzero provided that the past light-cones of summits x i and x j overlap at the initial time (or at least approach each other within distances m −1 ), as illustrated in the

Generalizations
Let us list here several extensions, for which the result of eq. (5.27) would remain valid modulo some minor changes: • Although we have assumed for simplicity in the derivation that all the observables are evaluated at the same time t f , the final result (5.27) remains valid for measurements at more general spacetime locations x 1 , · · · , x n . The only limitation is that all the separations between these points should be space-like, (x i − x j ) 2 < 0 if i = j, in order to avoid that the measurement performed at one point influences the outcome of the measurement performed at another point.
• It is possible to evaluate correlation functions where different observables are evaluated at each point x i , by having one type of node for each observable in the tree representation of eq. (5.27).
• One can easily replace the initial vacuum state by any coherent state. Instead of setting the initial field Φ ini (and its first time derivative) to zero after evaluating the derivatives corresponding to the links in the trees of eq. (5.27), one would have to set them to the values of Φ ini and ∂ 0 Φ ini that correspond to the coherent state of interest (see the appendix B for more details). A large initial value Φ ini of the coherent field only ensures that φ 1 . For the quasi-classical approximation to be valid to all orders in z(x), instabilities are necessary, like in the case of an initial vacuum state.
• Another extension is to consider a mixed state as initial state. If this state is highly occupied, the quasi-classical approximation is automatically satisfied to all orders in z(x), as explained in the appendix C.

Beyond the quasi-classical approximation
In the section 4, we have obtained the complete tree level result for the 1-point (eq. (4.9)) and 2-point (eq. (4.23)) functions, and one readily sees that they coincide with the result of JHEP09(2017)055 the quasi-classical approximation derived in the previous section (eq. (5.27) for n = 1 and n = 2, respectively). However, as we shall see now, for the 3-point function and beyond, the quasi-classical approximation does not include all the tree level contributions. This is expected from the fact that the quasi-classical approximation neglects some terms in the equations of motion for the fields φ 1 and φ 2 . In this section, we work out the full tree level result for the 3-point function, in order to clarify which terms are missed by the strong field approximation. Firstly, let us recall here the result for C {123} in the quasi-classical approximation: Let us now calculate in full the 3-point function at tree level, including the contributions that are beyond the quasi-classical approximation. To that effect, we must return to the original equations of motion (3.10), and expand them to second order in z, which leads to where we have systematically used the fact that φ 1 ≡ 0 in order to eliminate a few terms. The underlined term in the second equation is the only one that comes from the φ 3 1 φ 2 interaction term in eq. (3.5), that we had neglected in the quasi-classical approximation. The boundary conditions obeyed by these second-order coefficients at the final time read while those at the initial time relate their Fourier coefficients as follows Let us recall now that the quasi-classical approximation is exact for all the coefficients of lesser order, i.e. φ  1 does not contain any term coming from the vertex φ 2 φ 3 1 , its solution is identical to the result of the quasi-classical approximation. Let us now focus on the equation for φ (2) 2 . It has the structure of a linear equation of motion with the terms in the right hand side playing the role of source terms, since they do not contain φ (2) 2 itself. Therefore, we may decompose the solution as the sum of two terms; a term that solves the homogeneous (i.e. without source) equation and obeys JHEP09(2017)055 the non-trivial boundary conditions (6.4), and a term that solves the full equation with a trivial null initial condition: φ and Concerning ψ 2 , the equation of motion does not contain any term coming from the φ 2 φ 3 1 vertex, and all the objects that appear in the equation of motion and boundary conditions are exact in the quasi-classical approximation. Therefore, ψ 2 is itself correctly given by this approximation. For ξ (2) 2 , we can write the following solution: 1 (y; x 2 ) , (6.8) where G R (x, y) is the retarded propagator dressed by the background field φ 2 . The first term is already contained in the quasi-classical approximation since it does not involve the φ 2 φ 3 1 vertex, but the second term is present only beyond this approximation. Therefore, the only tree level term in the 3-point correlation function that is not contained in the strong field approximation is where we have used the explicit form (4.14) of φ 1 in order to rewrite it in a completely symmetric form in the second equality. The peculiarity of this contribution is that the correlation is created in the bulk by the self-interactions of the fields, instead of the pairwise initial state correlations that we have encountered in the quasi-classical approximation.  Generalizing the diagrammatic representation of eq. (5.27), such a term may be represented as follows: (6.10) (Note that in this representation, the power of the classical field φ (0) 2 (y) that accompanies the vertex does not appear explicitly.) The causal structure of this term is illustrated in the figure 5: because of the three retarded propagators, the vertex y must be located inside the overlap of the three past light-cones of summits x 1,2,3 , but it is not tied to the initial surface.
To close this section, let us compare the orders of magnitude of the contributions of eq. (6.1) and of eq. (6.10) in the presence of the instabilities discussed in the section 5.3. On the one hand, we have where the factor g 4 is relative to the completely disconnected 3-point function (each derivative δ/δΦ ini brings a factor g when the background field Φ is of order g −1 ). Regarding the contribution of eq. (6.10), each link brings a factor g since it replaces a field Φ, and the factor g 2 φ (0) 2 is of order g. Regarding the time dependence, each retarded propagator behaves as Therefore, we have We see that the terms beyond the quasi-classical approximation have the same order in g, but are exponentially suppressed by a factor of order exp(−µ(t f − t i )).

JHEP09(2017)055 7 Summary and conclusions
In this paper, we have studied the correlation function between an arbitrary number of observables measured at equal times (or more generally at points with space-like separations) in quantum field theory. Our main focus has been to obtain expressions that are valid in the presence of strong fields, which happens for instance when the fields are driven by a large external source. Firstly, we have constructed a generating functional that encapsulates all these correlation functions by coupling the observables to a fictitious source z(x). Using the retardedadvanced basis of the in-in formalism, we have shown that it can be expressed at tree level in terms of a pair of fields that obey coupled equations of motion and non-trivial boundary conditions that depend on the observables. At the first two orders in the fictitious source (i.e. for the 1-point and 2-point correlation functions), the results can be expressed in terms of the retarded classical solution of the field equation of motion, and its functional derivative with respect to its initial condition. Physically, this means that the 2-point correlation function at tree level is entirely inherited from the initial state fluctuations.
Beyond the 2-point function, the expansion in powers of z(x) becomes very cumbersome. However, in situations where the classical field dynamics has instabilities, a stronger quasi-classical approximation becomes legitimate, that leads to simplified field equations of motion that one may solve exactly. By doing so, one obtains an implicit equation for the generating functional of the correlation functions. From this functional relation, we have found that all the correlation functions are expressible in terms of the classical field and its derivatives with respect to the initial condition, generalizing the results for the 1-point and 2-point functions. Moreover, the expressions can be systematically represented by trees, where the nodes are the observables and the links are pairs of derivatives with respect to the initial condition of the classical field. Physically, these links correspond to correlations induced by fluctuations in the initial state. The reason why instabilities lead to this simple structure for correlation functions, in which correlations generated by interactions throughout the evolution are absent, is that the unstable modes generate a relative amplification (exponentially growing with time) of the initial state fluctuations versus the fluctuations in the bulk.
In the last section, we departed from the quasi-classical approximation and considered the tree-level 3-point correlation function in full. We found an additional tree-level term that does not exist in the quasi-classical regime, corresponding to correlations created in the bulk by the interactions themselves. We expect that such terms (and more complicated ones) exist in all n-point tree-level correlation functions for any n ≥ 3, beyond the quasiclassical regime.

JHEP09(2017)055
A Some properties of the mode functions In the section 4.3.2, we have introduced a basis of solutions for the linear space of solutions for a partial differential equation of the form where ϕ(x) is some real background field. This basis was defined as a set of solutions {a ±k (x)} labeled by a momentum k, whose initial condition is a plane wave of positive (in the case of a +k ) or negative (in the case of a −k ) frequency: Any solution of the equation (A.1) is a linear superposition of the functions a ±k . Given a solution a(x), let us define the following two components vectors Then, from two solutions a 1 and a 2 of eq. (A.1), we may define the following inner product reminiscent of the Wronskian for solutions of ordinary differential equations. This product is Hermitian, a 2 a 1 = a 1 a 2 * , (A.5) and constant in time: Its value is therefore completely determined by the initial conditions for the solutions a 1 and a 2 . In the case of the solutions a ±k introduced above, an explicit calculation gives a +k a +k = (2π) 3 2E k δ(k − k ) , A generic solution a of eq. (A.1) can be decomposed as follows a = d 3 k (2π) 3 2E k γ +k a +k + γ −k a −k , (A. 8) and from the inner products between the mode functions we readily see that the coefficients of this linear decomposition are given by γ +k = a +k a , γ −k = − a −k a . (A.9)

JHEP09(2017)055
Note that in the decomposition (A.8), the coefficients γ ±k are constant, and the time dependence is carried by the mode functions a ±k . Therefore, we can write a = d 3 k (2π) 3 2E k a +k a +k a − a −k a −k a . (A.10) Since this is true for any solution a, the following relationship must in fact be true Reinstating coordinates, this identity reads k a +k (x)ȧ −k (y) − a −k (x)ȧ +k (y) a −k (x)a +k (y) − a +k (x)a −k (y) a +k (x)ȧ −k (y) −ȧ −k (x)ȧ +k (y)ȧ −k (x)a +k (y) −ȧ +k (x)a −k (y) it is easy to check the following a in (p) χ = χ(p) χ , 3) The first equation tells us that χ is an eigenstate of annihilation operators, which is another definition of coherent states, and the second one provides the value of the normalization constant. Consider now the generating functional for the in-in formalism in this coherent state,

JHEP09(2017)055
where η(x) is a fictitious source that lives on the closed-time contour C introduced after eq. (2.5). The first step is to factor out the interactions as follows: .

(B.5)
A first use of the Baker-Campbell-Hausdorff formula enables one to remove the path ordering, giving where θ c (x 0 − y 0 ) generalizes the step function to the ordered contour C. Note that the factor on the second line is a commuting number and thus can be removed from the expectation value. A second application of the Baker-Campbell-Hausdorff formula allows to normal-order the first factor. Decomposing the in-field as follows, we obtain The factor of the first line can be evaluated by using the fact that the coherent state is an eigenstate of annihilation operators: . (B.9) We denote Φ χ (x) the field obtained by substituting the creation and annihilation operators of the in-field by χ * (k) and χ(k) respectively. Note that this is no longer an operator, but a (real valued) commuting field. Moreover, because it is a linear superposition of plane waves, this field is a free field:

JHEP09(2017)055
The second and third factors of eq. (B.8) are commuting numbers, provided we do not attempt to disassemble the commutators. Using the decomposition of the in-field in terms of creation and annihilation operators, and the canonical commutation relation of the latter, we obtain in (x), φ which is nothing but the usual bare path-ordered propagator G 0 c (x, y). Collecting all the factors, the generating functional for path-ordered Green's functions in the in-in formalism with an initial coherent state reads We see that it differs from the corresponding functional with the perturbative vacuum 8 as initial state only by the second factor, that we have underlined. This generating functional is also equal to 9 The first factor amounts to shifting the fields by Φ χ (x). The simplest way to see this is to write φ ≡ Φ χ + ζ . (B.14) In the definition (B.4), this leads to where the second factor in the right hand side is the generating functional for correlators of ζ. Comparing with eq. (B.13), we see that the generating functional for ζ is identical to the vacuum one, except that the argument φ of the interaction Lagrangian is replaced by Φ χ + ζ: