Tree-level correlations in the strong field regime

We consider the correlation function of an arbitrary number of local observables in quantum field theory. We show that, at tree level in the strong field regime, these correlations arise solely from fluctuations in the initial state. We obtain the general expression of these 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 from the strong field approximation, generated throughout the evolution of the system.


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 , 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 (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 possible, we consider a theory with a single real scalar field φ, whose Lagrangian density is given by 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 (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 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 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 Eq. (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) 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) 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 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 next-to-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 strong field 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 for 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 useful to discuss the strong field regime. 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 strong field 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 strong field 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.
2 Generating functional for local measurements

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) 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. 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: The operators are ordered 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. (13), 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. (13) 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 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)).
following explicit expressions: 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, these propagators 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) are given by where More explicitly, we have : (The vertices not listed explicitly here are obtained by permutations.) Finally, 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 should vanish when a measurement is performed.
3 First derivative at tree level 3.1 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. (26) 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 φ 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. (28)).

Classical equations of motion
Using the fact that G 0 12 and G 0 21 are Green's functions of + m 2 , respectively 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 at t f . For later reference, let us 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 (34) are easier to handle than the integral equations (28), 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 , (39) Note that the boundary term vanishes at the initial time t i , because G 0 12 is the retarded propagator. Likewise, we obtain the following equation for φ 2 : The boundary conditions at t i and t f are obtained by comparing eqs. (28) and (39-40). At the final time t f , the boundary condition is 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 −+ (see eqs. (14)), we see that, at the initial time, the combination φ 2 + 1 2 φ 1 has no positive frequency components, and the combination φ 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 (34), accompanied by the boundary conditions (41) and (45), 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 extremely hard 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: 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. (48), 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 φ 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 φ only contain the previously obtained φ (n ) 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 (34) and in the boundary conditions (41-43). This leads to considerable simplifications. Firstly, the boundary condition at t f for the coefficient φ From the equation of motion for φ Then, the boundary conditions at t i tell us that and since φ i.e. the usual Euler-Lagrange equation for the theory under consideration. Therefore, the coefficient φ 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 Since this object will appear repeatedly in the following, let us introduce a simple graphical representation: The shaded blob in fact encapsulates an infinite series of tree diagrams, a glimpse of which is given by the second equality.

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

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. (55). 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 φ can be expressed as a linear combination of these mode functions as follows, That this function satisfies the required equation of motion is obvious from the first of eqs. (58). The boundary conditions at the final time follow from the properties (122) of the mode functions that are derived in the appendix A (up to a factor i, the underlined combination of mode functions is the advanced propagator dressed by the background field φ (0) 2 ). The boundary conditions at the initial time tell us that the positive frequency content of φ 1 , and its negative frequency content is plus one-half that of φ (1) 1 . 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 φ 2 (t f , x 1 )) .
(60) Note that the underlined combination, that we denote G 22 , is nothing but the symmetric 22 propagator dressed by the background field φ (0) 2 . The 2-point correlation function is then given by 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. (61) 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. (61) can also be rewritten in a different way, that we will generalize to the case of the n-point correlation in the strong field approximation. Firstly, let us generalize the classical field φ 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 : 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 Thus, the tree level 2-point correlation function can be written as .
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 two derivatives with respect to the initial condition of the classical field Φ, one acting on each end of the link. We will represent diagrammatically this link 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 5 From its definition, the ⊗ product is symmetric, diagrams, made of two copies of the graphs that appear in eq. (54) 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 1-point 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 3point 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, this remains valid in the strong field approximation.
5 Correlations in the strong field regime 5.1 Strong field approximation Until now, our counting was based on the fact that a large external source J leads to large fields φ ± , but no 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 consider in this section an approximation that allows a formal solution to all orders in z. Here, we give only a very sketchy motivation for this approximation, and a lengthier discussion of its validity will be provided later in this section (after we have derived expressions for the fields φ 1 and φ 2 ).
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 interferences when squaring amplitudes, a quantum effect controlled by . Consequently, we may expect the difference φ 1 ≡ φ + − φ − to be small compared to φ ± themselves, i.e.
In this regime, that we will call the strong field approximation (SFA), we can approximate the equations of motion 6 (34) by keeping only the lowest order in φ 1 . This amounts to keeping only the terms linear in φ 1 in eq. (29) (in the case of a φ 4 theory, it means dropping the φ 3 1 φ 2 term in eq. (30)). In the approximation, they read while the boundary conditions are still given by (41) and (43). The problem one must now solve is illustrated in the figure 2. The field φ 1 obeys a linear equation Figure 2: Relationship between the fields φ 1 and φ 2 in the strong field approximation.
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. (71)) is formally identical to the first of eqs. (55). We can therefore mimic eq. (59) and write directly φ 1 as follows: where the mode functions a ±k should now be defined with φ 2 as the background, rather than φ 2 . To obtain eq. (72), it was crucial to have a linear equation of motion for φ 1 , a consequence of the strong field 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. Beside 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. (71). 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. (73) inside the exponential, this leads to Setting x 0 = t f and denoting 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. (73). 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. (78) 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 strong field approximation
Let us now return on the condition φ 1 φ 2 , that was used in the derivation of eq. (78), in order to see a posteriori when it is satisfied. To that effect, we can use eq. (72) for φ 1 . For φ 2 , the initial condition at t i is given by eq. (73).
For the sake of this discussion, it is sufficient to use a linearized solution for φ 2 in the bulk, that reads First of all, a comparison between eqs. (72) and (79) 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 when the classical solutions of the field equation of motion (52) 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. (79) is bilinear in the mode functions, we expect that Estimating the magnitude of φ 1 requires more care. Indeed, from eqs. (122) in the appendix A, antisymmetric combinations of the mode functions at equal times remain of order 1 even if individual mode functions grow exponentially with time. Thus, at the final time, we have 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. (72) is the advanced propagator G A in the background φ 2 . This retarded propagator may also be expressed in terms of a different set of mode functions b ±k defined to be plane waves at the final time t f , 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 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. (78) in powers of z
Although eq. (78) 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: and These examples generalize to all orders in z: the functional D[x 1 ; z] can be represented as the sum of all the rooted trees (the root being the node carrying the fixed point x 1 ) weighted by the corresponding symmetry factor 1/S(T ): (90) Note that the sum of the weights for the trees with n + 1 nodes (one of them being the root node x 1 ), is the n-th Taylor coefficient of the function w(z), that satisfies the following identity This functional identity may be viewed as a structureless version of eq. (78), in which the function z(x) is replaced by a scalar variable z. Eq. (93) 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 trees with n distinct labeled nodes 7 (including the node at x 1 ). Moreover, since derivatives commute, these successive differentiations eliminate the symmetry factors, leading to The number of trees contributing to this sum is equal to n n−2 . Eq. (95) tells us that, at tree level in the strong field 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 (95) 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. (96) that, after some hefty combinatorics, one may rewrite in the compact form of the left hand side of eq. (96).
In the strong field 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. (95), 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 figure 4.

Generalizations
Let us list here several extensions, for which the result of eq. (95) 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 (95) 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. (95).
• 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. (95), 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).
• Another extension is to consider a mixed state as initial state. If this state is highly occupied, the strong field approximation is also satisfied, as explained in the appendix C.

Beyond the strong field approximation
In the section 4, we have obtained the complete tree level result for the 1point (eq. (54)) and 2-point (eq. (68)) functions, and one readily sees that they coincide with the result of the strong field approximation derived in the previous section (eq. (95) for n = 1 and n = 2, respectively). However, as we shall see now, for the 3-point function and beyond, the strong field approximation does not include all the tree level contributions. This is expected from the fact that the strong field 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 strong field approximation: Let us now calculate in full the 3-point function at tree level, including the contributions that are beyond the strong field approximation. To that effect, we must return to the original equations of motion (35), 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. (30), that we had neglected in the strong field approximation. The boundary conditions obeyed by these secondorder coefficients at the final time read while those at the initial time relate their Fourier coefficients as follows Let us recall now that the strong field approximation is exact for all the coefficients of lesser order, i.e. φ does not contain any term coming from the vertex φ 2 φ 3 1 , its solution is identical to the result of the strong field 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 the non-trivial boundary conditions (101), 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 strong field approximation. Therefore, ψ (2) 2 is itself correctly given by this approximation. For ξ (2) 2 , we can write the following solution: where G R (x, y) is the retarded propagator dressed by the background field φ 2 . The first term is already contained in the strong field 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 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 strong field approximation. Generalizing the diagrammatic representation of eq. (95), such a term may be represented as follows: (Note that in this representation, the power of the classical field φ : Causal structure of the extra contribution to the 3-point correlation function that arises beyond the strong field approximation. In this contribution, the correlation is produced in the bulk, from the vertex φ 3 1 φ 2 . Causality implies that this vertex must be located in the intersection of the three light-cones. The red dot contains a power of the classical field φ 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. (98) and of eq. (107) 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. (107), 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 strong field approximation have the same order in g, but are exponentially suppressed by a factor of order exp(−µ(t f − t i )).

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 the strong field regime, that arises for instance when the fields are driven by a large external source and their classical equations of motion are subject to instabilities. 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 retarded-advanced 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.
For these two lowest orders, the result we have obtained is in fact exact at tree level, and does not rely on having strong fields. Beyond these low orders, this direct approach is very cumbersome to extend systematically. In order to circumvent this difficulty, we have introduced the strong field approximation, thanks to which one can formally solve the equations of motion with the appropriate boundary conditions, in the form of an implicit functional equation for the first derivative of the generating functional. 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.
In the last section, we departed from the strong field 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 strong field 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 strong field regime.

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 (111) 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. (111), we may define the following inner product reminiscent of the Wronskian for solutions of ordinary differential equations. This product is Hermitian, 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 generic solution a of eq. (111) can be decomposed as follows and from the inner products between the mode functions we readily see that the coefficients of this linear decomposition are given by Note that in the decomposition (118), the coefficients γ ±k are constant, and the time dependence is carried by the mode functions a ±k . Therefore, we can write Since this is true for any solution a, the following relationship must in fact be true Reinstating coordinates, this identity reads These identities are justification for eqs. (59) and (72).

B In-in formalism for an initial coherent state
A coherent state can be defined from the perturbative in-vacuum as follows where χ(k) is a function of 3-momentum and N χ a normalization constant adjusted so that χ χ = 1. From the canonical commutation relation it is easy to check the following a in (p) χ = χ(p) χ , 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, where η(x) is a fictitious source that lives on the closed-time contour C introduced after eq. (13). The first step is to factor out the interactions as follows: . (127) 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: 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: The second and third factors of eq. (130) 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 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 φ ≡ Φ χ + ζ .
In the definition (126), this leads to where the second factor in the right hand side is the generating functional for correlators of ζ. Comparing with eq. (135), we see that the generating functional for ζ is identical to the vacuum one, except that the argument φ of the interaction Lagrangian is replaced by Φ χ + ζ: In other words, the field ζ appears to be coupled to an external source and to a background field. Note that in the in-in formalism, the change φ → Φ χ + ζ applies equally to the fields φ ± on both branches of the time contour. Therefore, for the fields φ 1,2 in the retarded-advanced basis, we have and the integral equations (28) that determine these fields at tree level become Using φ 2 = Φ χ + ζ 2 and the fact that Φ χ is a free field, we see that the equations of motion corresponding to these integral equations are the same as eqs. (34). The boundary conditions at the final time read: while at the initial time we have the following relationship among the Fourier coefficients. The derivation of the main result (78) for a vacuum initial state can be reproduced almost identically for a coherent initial state, leading to the only difference being that the initial value of the classical field Φ is set to Φ χ instead of zero after the derivatives with respect to Φ ini have been performed.
C In-in formalism for a Gaussian mixed state Let us consider in this section a mixed initial state described by the following density matrix In this definition, β k has the same function as an inverse temperature (by allowing it to be momentum dependent we can also consider out-of-equilibrium systems). The expectation value in the pure initial state of eq. (9) is replaced by a trace 0 in · · · 0 in → Tr ρ · · · Tr ρ .
How to handle this type of initial state is well known from quantum field theory at finite temperature. The perturbative rules are identical to those exposed in the figure 1, but the propagators of eqs. (14) should be replaced by where we have defined In the retarded-advanced basis, the propagators G 0 12 and G 0 21 are unmodified, but the propagator G 0 22 becomes f k -dependent. The integral equations (28) are unaltered (but they acquire a hidden dependence on f k through G 0 22 ), and consequently the equations of motion (34) are unchanged. The boundary conditions (41) at the final time remain the same, but those at the initial time (45) are modified into φ (+) In other words, the factor 1 2 of eqs. (45), that can be interpreted as the zero point occupation of the vacuum, is now replaced by the total occupation number 1 2 + f k of the mixed initial state under consideration. We see from eqs. (148) that a large occupation number enhances φ 2 with respect to φ 1 . This is another situation where the strong field approximation introduced in the section 5 is applicable.