Loop diagrams in the kinetic theory of waves

Recent work has given a systematic way for studying the kinetics of classical weakly interacting waves beyond leading order, having analogies with renormalization in quantum field theory. An important context is weak wave turbulence, occurring for waves which are small in magnitude and weakly interacting, such as those on the surface of the ocean. Here we continue the work of perturbatively computing correlation functions and the kinetic equation in this far-from-equilibrium state. In particular, we obtain the two-loop kinetic equation for waves with a cubic interaction. Our main result is a simple graphical prescription for the terms in the kinetic equation, at any order in the nonlinearity.


Introduction
Kinetic equations are widely used to study equilibrium and non-equilibrium dynamics, transport, and even turbulence.The kinetic equations are a truncation of a vast hierarchy, involving multimode correlators.An important problem is understanding when this truncation is correct.It has long been recognized in the context of particles (the Boltzmann equation) that corrections can be large and give qualitatively new effects [1][2][3][4], and have been studied in a number of works [5][6][7][8][9][10][11].
The kinetics of waves [12] is particularly fascinating as it allows for the study of weak wave turbulence [13], an old subject [14] undergoing a recent experimentally driven revival in a range of physical contexts, see e.g.[15] and references therein, going beyond the canonical example of ocean waves [16,17].For waves, unlike for particles, there hasn't been a systematic study of higher order corrections to the leading order kinetic equations.
The connection between kinetic equations and quantum field theory [18] has given an effective way of computing corrections, and made it evident that they are important, modifying the standard assumptions on the range of applicability of wave turbulence [19].In short, even if the coefficients of the nonlinear terms in the Hamiltonian are numerically small, intermediate states in the scattering of waves (loop diagrams) can have interacting waves with vastly different momenta.This large ratio of momenta can then overpower the smallness of the interaction coefficient, making higher order loop diagrams seemingly dominant rather than suppressed.Summing an infinite class of loop diagrams may then give a new, renormalized, kinetic equation, which encodes phenomenon not captured by the standard leading order kinetic equation [19].The need to compute loop diagrams, in a simple and compact way, and potentially to high orders, is acute.This is what we seek to address, extending our previous work [18,20].Our discussion will be confined to classical kinetic theory (although we study it using tools of quantum field theory); the extension to quantum kinetic theory will be discussed in [21] and is relevant for studies of thermalization and wave turbulence in QCD, see e.g.[22][23][24][25][26].
In Sec. 2 we review the context of interacting waves and the leading order in nonlinearity (tree-level) correlation functions and the corresponding kinetic equation.In Sec. 3 we compute the one-loop correction to the correlation functions.We do this for a theory with a cubic or quartic interaction.In Sec.3.3 we use this to give the next-to-leading order kinetic equation for a theory with a cubic interaction.This extends the results of [18,20] which studied the simplest case, of a quartic interaction.In Sec. 4 we give our main result: a prescription for almost immediately writing down the contribution of any Feynman diagram, at any order in the coupling, to an equal time correlation function.This turns the problem of finding higher order terms in the kinetic equation into a simple bookkeeping exercise.We conclude in Sec. 5.
In Appendix A we show that an alternate method for deriving the higher order terms in the kinetic equation, by perturbatively solving the Liouville equation and performing a phase space average [20], exactly reproduces -at each order in perturbation theory -the method used in the main body, of averaging over external Gaussian random forcing which is set to zero at the end.In Appendix B we discuss symmetries of the Hamiltonian and how it relates different terms in the kinetic equation, which is discussed further in Appendix C. A few technical remarks are relegated to Appendix D.

Interacting waves
Our context will be that of classical interacting waves.This occurs in many situations, ranging from surface gravity waves to spin waves.We focus on a single field, which can, for instance, describe the height of the ocean.The variables are the field and its canonical conjugate, though it is often simpler to work with a single complex field a p , which is the sum of these two.The most general Hamiltonian is a sum of q-body Hamiltonians, (2.1) In particular, H 2 is the free Hamiltonian, where a † p denotes the complex conjugate of a p , a † p ≡ a * p .Generally, if one is considering a theory with weak interactions, the higher order the interaction term the more suppressed it is.For this reason, it is common to focus on solely the lowest order nonvanishing term.In general, the first interaction term that dominates is the cubic term, so that the Hamiltonian is truncated to, where for the interaction term there is an implicit momentum conserving delta function δ(⃗ p 1 −⃗ p 2 −⃗ p 3 ).
In sections that follow, it will be convenient to use the shorthand λ 123 ≡ λ p 1 p 2 p 3 , and similarly for larger q.In some cases the dispersion relation is such that a resonant cubic interaction is forbidden, i.e., it is not possible for the wave to have both ⃗ p 1 = ⃗ p 2 + ⃗ p 3 and ω p 1 = ω p 2 + ω p 3 .Then the lowest interaction term is quartic.This is the case we focused on in [18] and [20] (the symmetries of the quartic term make it slightly simpler than the cubic term), 1 In the main body of the text we will focus on the case of q = 3 and q = 4.The generalization to arbitrary q is straightforward.

Averaging over forcing
Since this is a many-body chaotic system, one has to perform some kind of average in order to have sensible quantities.The averaging that we do here is to add Gaussian-random forcing.In particular, the equations of motion, are modified so as to add forcing, as well as dissipation in order to absorb the flux of energy.The new equations of motion are 2 , 1 The cubic couplings have the symmetry λ p 1 p 2 p 3 = λ p 1 p 3 p 2 and the quartic couplings have the symmetry There are other choices of dissipation one could make, such as In the limit of vanishing dissipation these are all equivalent.
where the forcing is drawn from a Gaussian distribution with variance F k , After computing the correlation functions, we take F k , γ k → 0 with fixed n k , The reason for this notation is that n k is the occupation number of mode k for the noninteracting theory.Going forward, for notational simplicity, instead of writing F k , we will write 2γ k n k .
It is important to note that we will be computing correlation functions in a stationary state.
In other words, we assume that we have a stationary state and then later self-consistently pick F k and γ k to ensure this.It will turn out that there are multiple possible stationary states in the zero dissipation limit: the thermal state and the turbulent state.However, nothing in our calculations of the correlation functions and kinetic equation is dependent on the properties of these states; we only need to assume we are computing correlators about a stationary state.
One can compute correlation functions by solving the equations of motion perturbatively in the interaction and then averaging over forcing.For each term in the perturbative expansion one can associate a corresponding Wyld diagram [28,29].
A more streamlined method for doing this calculation was introduced in [18], and consists of integrating out the forcing at the outset.In particular, our classical stochastic field theory is equivalent to a quantum field theory, with expectation values given by a path integral for a Lagrangian that is the square of the classical, force-free, equations of motion, Let us break up the Hamiltonian into the free part, H 2 , and the interacting part, where where c.c. denotes the complex conjugate.The term contained in the dots of L int that we have left out is the square of the absolute value of δH int δa k ; we can neglect this term provided that when evaluating Feynman diagrams we drop contact interactions (when two times collide), see Appendix

D.
Let us work out the Feynman rules.As usual, the Feynman rules are most convenient in momentum/frequency space.The quadratic term in the Lagrangian, L f ree , gives the propagator, (2.12) We will use shorthand D i ≡ D p i ,ω i , n i ≡ n p i , and γ i ≡ γ p i .Notice that in the limit of γ k → 0, the propagator becomes a delta function, n k 2πδ(ω − ω k ).It will, however, be important to keep γ k finite until the very end of the calculation.It will be convenient to define,

.13)
The Feynman rule for the vertex is3 where we have introduced the notation ω 12;3 ≡ ω 1 +ω 2 −ω 3 and ω 12;34 ≡ ω 1 +ω 2 −ω 3 −ω 4 .The treelevel frequency space correlation functions are trivially obtained by adding external propagators and accounting for the appropriate combinatorial factors from Wick contractions, The corresponding Feynman diagrams are shown in Fig. 1.The time-space correlation functions are Fourier transforms, In particular, the equal-time tree-level correlation functions are, ) where ϵ = i γ i represents the sum over the dissipation constants appearing in the vertex.
An important quantity is the occupation number n k of mode k, n k = ⟨a † k a k ⟩, which is governed by the kinetic equation.Through use of the equations of motion it can be expressed in terms of the equal-time correlation function.In particular, for q = 3, where all operators on the left and the right are evaluated at time t.For the general Hamiltonian (2.1) this has the obvious generalization.
Inserting the tree-level correlation function (2.18) into (2.20)gives the standard kinetic equation for waves [13].Our goal will be to compute higher order in the coupling corrections to this.In the next section we will therefore turn to computing the equal-time correlation functions perturbatively in the coupling.

Loops
The tree-level correlation functions were given in the previous section.Here we will compute their one-loop corrections, first in the case of the quartic interaction and then for the cubic interaction.

Quartic interaction
The one-loop (order λ 2 ) correction to the tree-level four-point function for the case of q = 4 (quartic interaction) was computed in [18].One of the corresponding diagrams is shown in Fig. 2; there are two more diagrams with arrows in different directions which we won't discuss here.Here we redo the calculation in a faster way, working in time space instead of frequency space, which is better suited for generalization to higher order loops.
Using the Feynman rules, the contribution of this diagram to the one-loop four-point function,  We exchange the order of integration, and first evaluate the ω i integrals.The ω i integrals are trivial to evaluate, by closing the contour and picking up simple poles.For any particular choice of time orderings of t a , t b , t c , convergence imposes a unique choice of if the contours of integration should be closed in the upper or lower half-plane.In particular, the poles come from the propagators D i in V .For the integrals dω 2π e −iωt ij : if t ij > 0 then we close in the lower-half plane, picking up the pole ω = ω p −iγ.If t ij < 0 then we close in the upper half plane, picking up the pole ω = ω p +iγ. 5We split the time integrations into 3! regions, corresponding to the 3! possible time orderings of t a , t b , t c .Unless t a is the earliest time, one finds a vanishing contribution.This leaves two possible orderings: We close the ω 1 , ω 2 integrals in the lower half plane, and close the ω 3 , ω 4 , ω 5 , ω 6 integrals in the upper half plane.The vertex factor V becomes, and the time integral becomes, dt b e −i(ω p 1 +ω p 2 −2iγ)t ba e i(ω p 3 +ω p 4 +2iγ)t ca e i(ω We close the ω 1 , ω 2 , ω 5 , ω 6 integrals in the lower half plane, and we close the ω 3 , ω 4 integrals in the upper half plane.This vertex factor becomes and the time integral becomes, dt c e −i(ω p 1 +ω p 2 −2iγ)t ba e i(ω p 3 +ω p 4 +2iγ)t ca e i(ω The sum of (3.5) and (3.8) reproduces what we found earlier in [18] (see Eq. 4.25 and use the identity A.8 in [20]).6

Cubic interaction
Consider now the equal-time three-point function ⟨a 1 a † 2 a † 3 ⟩(t a ) arising from a cubic interaction.The tree-level correlation function, at order λ, was given earlier in (2.18).At order λ 2 all Feynman diagrams vanish.At order λ 3 one has two kinds of diagrams: the "tetrahedron" diagram, which can be viewed as a renormalization of the cubic interaction vertex, and a loop diagram arising from renormalization of the propagator.
where the six terms arise from the 3! possible orderings of the times t b , t c , t d > t a .In particular, the time orderings of the six terms are: from first to last, respectively.In fact, as a result of the symmetries of the tetrahedron, each of the six terms can be obtained by a symmetry transformation of one another, see Appendix B.
In addition to this diagram, there is the same tetrahedron diagram but with the arrow on line 5 reversed.The result for the latter can be obtained from the former in a simple way, from symmetry, and is also discussed in Appendix B.

Propagator renormalization
Now we look at the loop diagram shown in Fig. 4 which can be viewed as arising from propagator renormalization.We find that its contribution to the equal-time three-point function is, where the terms correspond to the time orderings to last, respectively.7 There are several more diagrams of this kind, which can be viewed as a symmetry transformation of this diagram, as we will see below.In fact, in (3.10) one should, by momentum conservation, set p 4 = p 1 .This gives some terms that are divergent.Specifically, (3.10) gives the term, However, the term in brackets vanishes by the tree-level kinetic equation.This is a useful illustration of the fact that all our computations are of equal-time correlation functions, assuming a stationary state.

Next-to-leading order kinetic equation for cubic interactions
We are now ready to write down the kinetic equation to next-to-leading order.It follows from the equal-time three-point function, via (2.20).The complex conjugate version of this equation is, The three-point function consists of an infinite sequence of terms, grouped by their order in λ, and we go up to order λ 3 , The leading order term ⟨a p 1 a † p 2 a † p 3 ⟩ λ , given earlier in (2.18), gives the standard (leading order) kinetic equation.The next-to-leading order term ⟨a p 1 a † p 2 a † p 3 ⟩ λ 3 gives the first corrections to this.The terms appearing here are: V a , which arrises from the diagram shown earlier in Fig. 3(a We get P c and P d from P a and P b , respectively, by sending 1 ↔ −3 and 6 → −6 (and also 4 → −4, however 4 in a is just 1).Finally, we get diagram e from diagram a by sending 1 ↔ −3 and 6 → −6.Note that diagrams c, d, e have a propagator renormalization for line 3.We could have alternatively put it on line 2; we account for this by adding a factor of 2 for P c , P d , P 3 .
which get the propagator renormalization.Explicitly, Further discussion of terms in the kinetic equation is given in Appendix C.

Prescription for a general Feynman diagram
In this section we give an algorithm for computing the contribution of any Feynman diagram to an equal-time correlation function.In Sec. 2 we already showed that correlation functions can be computed order by order in the nonlinear interaction (λ), just as one does in quantum field theory.The remaining task is, for a given Feynman diagram, to perform the integration over the intermediate times (or, equivalently, frequencies) appearing in the loops.We did this explicitly in the previous section, for some particular one loop diagrams.From the several diagrams that we evaluated, we can deduce the rule for a general diagram.We first state the rules and then derive them.

Rules:
1. Pick an ordering of the times at each vertex.The time t a at which the correlation function is being evaluated must be the smallest time.
2. Start at the latest time on the diagram and move from vertex to vertex in decreasing order of their times, until finally reaching the vertex at the earliest time.Each next vertex must be a neighbor of at least one previously visited vertex.
3. At each step in this process, draw an imaginary loop enclosing all vertices visited so far.
Write down a factor of −1 where ω i , ω j , . . .are the frequencies of all lines entering this imaginary loop and ω a , ω b , . . .are the frequencies of all lines leaving the imaginary loop.
4. In addition, at each step in the process, write down a factor of where ω k , ω l , . . .are the frequencies of the lines entering the vertex and ω α , ω β are the frequencies of the lines leaving the vertex.However, omit an 1 n i if it has already been included at a previous step.Justification: Let us understand these rules.Suppose for concreteness that we have four vertices with the time ordering t d >t c >t b >t a .The integrals we need to do are, Here the σ are the sum of the outgoing minus ingoing frequencies from a vertex; we get these through the delta functions we insert, see below (3.1).One now sees that as we successively do the integrals, starting from the latest time, we get the rules listed in Steps 2 and 3. Note that all our integrals here converge (because the relevant poles of the ω integrals were chosen to ensure this), so the contributions from infinite time vanish.Now let us understand the vertex factor, appearing in Step 4. According to the Feynman rules, each vertex comes with a sum of g * i for each propagator leaving the vertex and a −g i for each propagator entering the vertex, see (2.14).At the poles, ω i = ω p i ± γ i , we see from (2.13) that g i and g * i become, So from each vertex we will get a sum and difference of various 1/n i .We need only determine if for a propagator entering or leaving a vertex we get a 1/n i or nothing.Consider two vertices next to each other, as in Fig. 2(a).In the corresponding vertex factor (3.1), either the g 5 +g 6 term coming from vertex b gives 1/n 5 +1/n 6 and the g * 5 +g * 6 term coming from vertex c gives 0, or vice-versa.Which option is the correct one depends on the time ordering of t b and t c , which determines if the ω 5 , ω 6 integrals are closed in the upper or lower half plane.In particular, the vertex at the later time gets to have the 1/n 5 +1/n 6 .Since in our prescription we proceed from later to earlier time, this explains Step 4. and Step 4 instructs us to write a factor of The next largest time is t c , so we expand our loop (shown in green) to include vertex c.
Step 3 instructs us to write a factor of −1 and Step 4 instructs us to write a factor of As Step 4 says, even though ω 1 and ω 2 are leaving the green circle, we do not include since this term already appeared earlier, when accounting for the blue loop around vertex b, (4.6). Implementing Step 5, we reproduce (3.5).The other time ordering, t c >t b >t a , proceeds in a similar fashion.We have drawn the corresponding loops in Fig. 6(b).
Let us now look again at the loop digram for the cubic interaction, Fig. 3 where there is no 1/n 2 −1/n 4 , since it was already included in (4.9).Finally, we expand the loop to the loop shown in red in Fig. 6(c), so as to include vertex c.Frequencies ω 2 and ω 3 are entering and ω 1 is leaving.Steps 3 and 4 instruct us to write a factor of where there is no 1/n 2 or 1/n 3 , since both were already included.Multiplying (4.9), (4.10), and (4.11), we reproduce the second of the six terms in (3.9).
As our final example, let us look at an example of a diagram that occurs if we have both cubic Applying the rules we get the contribution, where reading from left to right, the terms correspond to larger and larger loops shown in Fig. 7 (b).

Discussion
In [18] we gave a systematic method for computing higher order terms in the wave kinetic equation: correlation functions are defined by averaging over Gaussian random forcing, which is set to zero at the end.This stochastic classical field theory is equivalent to a quantum field theory.
One computes equal-time correlation functions in the quantum field theory, order by order in the number of interaction vertices, and inserts into the right hand side of (2.20) to get the kinetic equation.The only nontrivial step is evaluating the Feynman diagrams with loops, which involve integrals over the intermediate times (or frequencies) inside the loops.In this paper we have shown how to carry out the integrals in general, obtaining a simple prescription for writing the result, described in Sec. 4.
Our answer still has integrals over the loop momenta; these need to be done on a case by case basis, as they are functions of the couplings, λ p 1 p 2 p 3 and λ p 1 p 2 p 3 p 4 , which are momentum dependent and system dependent.If these loop integrals over momenta didn't have UV or IR divergences, we would be close to being done -the corrections to the leading order kinetic equation would always be small, and the largest corrections would come from the loop diagrams with the fewest number of interaction vertices.But we're not done: the loops often have UV and IR divergences [19].
This introduces additional small and large parameters which compete with the smallness of the nonlinearity parameter, and make the perturbation theory multiscale.We must decide which classes of loop diagrams are dominant and sum those.In other words, we must renormalize.The results in this paper provide the necessary tools to carry out this task, which will be the subject of future work.

A. Perturbative Liouville equation for waves
An alternate method for studying wave turbulence is to average over phase space, instead of introducing external forcing and dissipation.In particular, one aims to find the evolution of the phase space density ρ(J i , α i , t), a function of the action, J i , and angle, α i , variables, with i running over the different modes of the field.This method was used in [20] to find the next-to-leading order correction to the kinetic equation for a theory with four-wave interactions.Here we extend this to the theory with three-wave interactions, like the one studied in the main body, reproducing the results found there.In fact, we go further: we show that at every order in perturbation theory, the results of this method agree with the one in the main body of the text.8

Perturbative solution of the Liouville equation
The evolution of the phase space density is governed by the Liouville equation [11], where L 0 is due to the free part of the Hamiltonian, H 2 (which only depends on the action variables), and δL is due to the interacting part, H int .We have defined, ω j ≡ ∂H 2 ∂J j .The solution for a time-independent phase space density is given by iterating, where the eigenfunctions |⃗ n⟩ of the free Liouville operator L 0 are plane waves ⟨⃗ α|⃗ n⟩ = exp (i⃗ n • ⃗ α) and G(⃗ n) is the propagator, The kinetic equation is then given by, Inserting (A.3) gives δL = (δL) first + (δL) second + . .., where where ρ(J) ≡ ⟨0|Φ 0 ⟩.We will take ρ(J) to be a Gaussian, The main point is this: the phase space integrated over all angles, ρ(J), which has no angular dependence, is governed by a differential equation, the right-hand side of which encodes a series of transitions to intermediate states which have angular dependence.In particular, each interaction vertex gives q of the different modes angular dependence (for a q-body interaction).We start and end in a state with no angular dependence.
In [20] we specialized these equations to interacting waves with a quartic interaction.Here we look at these equations for a field theory with a cubic interaction.The Hamiltonian, H 2 + H 3 , takes the following form when written in action-angle variables, a p = J p e −iα p , Computing the Liouville operator we have that L 0 is (A.2) with ω i given by ω p and δL in (A.2) is equal to, where ∂ i ≡ ∂ J i and where ⃗ e i;j,k , which is a vector that has a −1 in the i'th entry and a +1 in the j'th and k'th, and a zero elsewhere; in other words, −⃗ e i;j,k • ⃗ α = α i −α j −α k .The matrix element of δL trivially follows, where ⃗ e 2,3;1 = −⃗ e 1;2,3 .

Leading order
We now start iteratively computing the kinetic equation (A.5), starting with (δL) first in (A.6).
There are two diagrams, as shown in Fig. 8.The first, Fig. 8 where we made use of the commutator ∂, an extra minus sign for ω p 1 ;p 2 ,p 3 .Adding the two and using (A.5) we recover the leading order kinetic equation,

Next-to-leading order
We now look at the terms that appear at next order.All terms in (δL) second (A.7), with two intermediate states vanish, so we look at (δL) third , which has three intermediate states.One such Using ρ( ⃗ J) (A.8) and evaluating the integral of δL t ({a, b, c, d}) multiplied by J r , as prescribed by (A.5), gives, dJ J r δL t ({a, b, c, d}) = This matches the last term in (3.9).In more detail, in the main body of the text, the kinetic equation resulted by inserting contributions to the equal-time three-point function, such as (3.9), into (2.20).The contribution to the kinetic equation of the last term of (3.9) is what essentially matches (A.16).Notice that (2.20) has a factor of (δ kp 1 −δ kp 2 −δ kp 3 ) in front, which is just the factor of (δ 1r −δ 2r −δ 3r ) in (A.16) in slightly different notation.
There are other possible intermediate states.We can obtain them by changing the order in Of course, both these diagrams are simply redrawings of the tetrahedron diagram that appeared in Fig. 3.We need to consider all 24 permutations of the vertices.The six terms with vertex a first will give the six terms in (3.9).The six terms with vertex c first will give the same thing, but complex conjugated and with a relative minus sign.So their sum is just the imaginary part of the terms with vertex a first.We therefore reproduce the contribution of (3.9) to the kinetic equation (notice that (2.20) takes the imaginary part of the coupling multiplying the three-point function).The terms with vertex b first combined with the terms with vertex d first (which are the complex conjugate), reproduce the contribution of the tetrahedron in Fig. 3) with 5 → −5.We may do a similar analysis for the propagator renormalization diagram, but this is unnecessary, as we will now give a general argument that the results will always match what is in the main body of the text.

General Prescription
In implementing the perturbative solution of the Liouville equation, (A.vertex, separating the edges entering it, then we have a diagram that contributes to a q-point correlation (for a q-body interaction).In particular, it may be a diagram for either a correlator such as ⟨a 1 a † 2 a † 2 ⟩ (for q = 3), or its complex conjugate ⟨a † 1 a 2 a 3 ⟩.The complex conjugate arrises if we reverse all the arrows on the Feynman diagram.We can account for this by including only one Feynman diagram out of each pair that have all arrows reversed relative to each other and then taking the imaginary part, as we saw earlier.We can achieve this by drawing all diagrams contributing to ⟨a 1 a † 2 a † 2 ⟩ (for q = 3) and connecting the external lines together to meet at a vertex.We then consider all possible orderings of the other vertices, which determines the intermediate states.
At this stage the number of terms we have matches what was discussed in Sec. 4. Now let us show that the terms are in fact the same.We get from vertex to vertex with a propagator , where |⃗ n⟩ is the current state.This matches the prescription given in Step 3 (4.1) of having a denominator with a sum of all incoming frequencies minus all outgoing frequencies.In our situation here the imaginary loop encloses all vertices visited so far.At each vertex we also have a transition matrix of δL, see (A.10) or (A.11).We see that it it involves where ω k , ω l , . . .are the frequencies of the lines entering the vertex and ω α , ω β , . . .are the frequencies of the lines leaving the vertex.The transition matrix element also has factors of √ J i and of n ′ i /J i , see (A.11).However, the latter can be eliminated by appropriately commuting through the The end result, as seen in, for example, (A.15), is that one has each J i appear once, immediately to the right of when the corresponding ∂ i appears.After inserting ρ(J) that is an exponential, (A.8), it is clear we reproduce Steps 4 and 5 (4.2).Similar rules were found in [30,31].

B. Symmetries
In this Appendix we show that the symmetries of the Hamiltonian, combined with the symmetries of the some of the Feynman diagrams, provide an efficient way of finding many of the contributions to the kinetic equation.We will, for instance, show that once one has any of the six terms appearing in the contribution of the tetrahedron diagram (3.9) the five others follow by a symmetry transformation.

Symmetries of the Hamiltonian
We first note the symmetries of the Hamiltonian, written in action-angle variables (A.9).The Hamiltonian is invariant under, Note that λ * 123 transforms in the same way, λ * 123 → −iλ 123 .The Hamiltonian is also invariant under a flip in sign of just two indices, or a flip in one index, This is useful because we sometimes have diagrams which have some arrows flipped relative to other diagrams.Flipping an arrow means flipping the sign of α.Due to the symmetries of the Hamiltonian we can instead flip the signs of the ω p i and J i , and change the couplings, as stated above.Note that in the kinetic equation we have n i instead of J i (the former is the expectation value of the latter).So an arrow going in the opposite direction on line i means that in the contribution to the kinetic equation n i , ω p i , and δ ir all pick up a minus sign.We mentioned in the main body that, in addition to the tetrahedron diagram in Fig. 3 that we evaluated, there is another diagram, which is just Fig. 3 with the arrow on line 5 reversed, as shown in Fig. 11.In our new notation, 5 → −5.This means that its contribution is (3.9) with n 5 → −n 5 , ω p 5 → −ω p 5 , and the couplings λ * 356 λ 145 replaced by −λ 635 λ * 451 .

Symmetries of the diagrams
We will now show that the different contributions of a single tetrahedron diagram are related by symmetry as well.

Transformations of a tetrahedron
Let us look at the tetrahedron diagram studied in Sec.3.2, and shown again below in Fig. 12(a), and consider all its symmetry transformation.There are 24 in total, corresponding to the 4! possible positionings of the vertices.The symmetry transformations that we use are rotations of the tetrahedron as well as mirror reflections.For instance, in Fig. 12 we have shown the 6 transformations that leave vertex a fixed: Fig. 12    This is recorded as the first entry in our table below.We see that, for instance, Fig. 12(e) has vertices a, b, d, c and edges −2, −1, 3, −4, 6, 5, where we are using a minus sign to denote the arrow running in the opposite direction.We have recorded this as the second entry in our table.The other entries are found in a similar manner.

vertex ordering edges
Transforming terms in the kinetic equation Let us now show how we can use these symmetries to relate different terms in the kinetic equation.In particular, we start with one of terms, (A.16), and define the right hand side of (A.16) without the couplings to be t(a, b, c, d), We view t(a, b, c, d) as a of the six edges, t(1, 2, 3, 4, 5, 6).The contribution to the kinetic equation of the other terms will still be given by t, but the arguments will contain some permutation of these six variables, and they may have minus signs.For instance, consider the term in which we have swap vertices c and d, which by our table corresponds to, We now collect all 6 terms that have vertex a first.We denote its contribution to the kinetic equation by, where T c can be related to T a by a rotation of the tetrahedron that transforms vertex a into vertex c.There are multiple ways of doing this, which differ by permutations of the other three vertices.
Picking the rotation {a, b, c, d} → {c, d, a, b}, and using the corresponding entry in the table, we There is no need to compute T b , because it can be related in a simple way to T a : we simply rotate vertex a into vertex b.There are multiple of doing this, which differ by permutations of the other three vertices.Picking the rotations {a, b, c, d} → {b, d, a, c}, and using the corresponding entry in the table, we have, Note the minus sign in front of (B.12) was placed there to account for minus sign acquired by the couplings under this transformation.Explicitly, T b is,

C. Manipulating the kinetic equation
In the main body of the text we found that the contribution of the tetrahedron diagram to the kinetic equation is, 2i Im where T a was given in (B.8) or equivalently the complex conjugate of (3.9).The second term is the tetrahedron diagram in Fig. 3 with the arrow on line 5 reversed, as discussed in Appendix B.
This form of the kinetic equation is acceptable, but the form of the answer that has a more clear physical interpretation is one in which there are explicit delta functions.In particular, we write, where the first term comes with an implicit principal value.Using this, we may group terms based on the number of delta functions.This is what we partially do in this appendix.
One consistency check for our kinetic equation is that the thermal state should be stationary.
In other words, inserting n i = 1/ω p i into (C.1)should give zero.This is not manifest for the piece of (C.1) that has no delta functions.We will rewrite (C.1) in a way that will make it manifest, by showing that there is no such piece.
Another way of writing (C.1) is as,

One delta function
Next, we look at the piece of C.1 which has one delta function.We will see that we can write this contribution to the kinetic equation entirely in terms of the following two functions, Like with T a , these are functions of (1, 2, 3, 4, 5, 6).There are a few special vertex orderings we will want to consider.We write them, along with the corresponding edges.The two functions we just wrote are viewed as being at the vertex, A ≡ {a, b, c, d} = (1, 2, 3, 4, 5, 6).We will also need, By the same manipulations as done in the beginning of this section, this is equal to 4 To summarize: we have taken the part of the kinetic equation due to the tetrahedron diagrams and extracted and simplified the piece that has one delta function.The result is the sum of (C.12) and (C.15).
One can continue, and extract the piece of the kinetic equation that has a product of two delta functions, and the piece that has a product of three delta functions.As we go to higher order in the nonlinearity, the kinetic equation will have products of an increasing number of delta functions.It would be useful to have a general prescription for extracting the coefficients of these delta functions.This would be a kind of extension of the cutting rules in finite temperature field theory, see e.g.[32].

D. Contact terms and lollipop diagrams
In this appendix, we provide justification for disregarding the square of the absolute value of δH int δa k in the interaction Lagrangian (2.11) (redundant interaction).Furthermore, we demonstrate that the existence of a non-zero expectation value for the complex field a p in the cubic theory (2.3) does not affect our findings.We begin by addressing the redundant interaction.The first term on the right-hand side is obtained by expanding the path integral to the mth order in the interaction displayed in (2.11).The ellipsis represent terms involving the square of the where the complex constant b 0 is defined by ⟨a p ⟩ = b 0 δ(⃗ p); it can be calculated by evaluating the diagrams in Fig. 13   where quantities with tilde absorb various terms which emerge due to the field redefinition (D.10). 9 By construction, the last two terms guarantee that ⟨a 0 ⟩ = 0.
We have shown that the expectation value of a p in the cubic theory can be ignored by appropriately redefining the frequencies, couplings, and fields.To achieve this, one introduces a linear term of the form J 0 a † 0 + J * 0 a 0 in the Hamiltonian (2.3), where the unspecified complex constant J 0 is appropriately tuned such that ⟨a p ⟩ = 0 to all orders of the perturbative expansion.In other words, we can safely disregard the lollipop diagrams in Fig. 13.The value of ⟨a p ⟩ is determined by minimizing the complete effective potential of the model, whereas the Hamiltonian (2.3) describes only excitations.

Figure 2 :
Figure 2: (a) A one-loop diagram contributing to the four-point function.(b) We set the external times to be equal to t a , represented by the four lines meeting at the point t a .

Figure 3 : 2 p 4 ,p 5 ,p 6 λ 426 λ * 356 λ 145 6 i=1
Figure 3: (a) A one-loop diagram contributing to the three-point function.(b) We set the external times to be equal to t a , represented by the three lines meeting at the point t a .We refer to this as the tetrahedron diagram

Figure 4 :
Figure 4: (a) A one-loop diagram contributing to the three-point function.(b) We set the external times to be equal to t a , represented by the three lines meeting at the point t a .

Figure 5 :
Figure5: The different propagator renormalization diagrams contributing to the three-point functions.All the diagrams can be obtained from a (whose contribution is denoted by P a ) through symmetry transformations, see Appendix B. In particular, we get P b from P a by sending 5 → −5.We get P c and P d from P a and P b , respectively, by sending 1 ↔ −3 and 6 → −6 (and also 4 → −4, however 4 in a is just 1).Finally, we get diagram e from diagram a by sending 1 ↔ −3 and 6 → −6.Note that diagrams c, d, e have a propagator renormalization for line 3.We could have alternatively put it on line 2; we account for this by adding a factor of 2 for P c , P d , P 3 .

5 .
Multiply the resulting expression by a product of n i , one for each frequency ω i appearing in the diagram, and multiply by a product of the couplings for each vertex, and sum the resulting expression over all internal momenta.6. Repeat Step 1 through 5 for all possible time orderings and sum the results.Include the appropriate Feynman diagram combinatorics factor.

Figure 6 :
Figure 6: Applying the rules to four-point and three-point correlation functions.

−1 ω p 2 ,p 6 ; 9 )
, which we evaluated earlier.There are six possible time orderings.Let us look at, for instance, the one in which t d >t b >t c >t a .We draw a loop, shown in green in Fig. 6(c), around vertex d, which is at the latest time.Frequencies ω 2 and ω 6 are entering and ω 4 is leaving.Steps 3 and 4 instruct us to write a factor of Expanding the loop to include the vertex b, which is at the next largest time, we get the loop shown in green in Fig. 6(c).Frequencies ω 2 and ω 3 are entering and ω 4 and ω 5 are leaving.Steps 3 and 4 instruct us to write a factor of

Figure 7 :
Figure 7: (a) An example of a higher order diagram.(b) Applying our rules to a particular time ordering, resulting in the contribution (4.12).
3), we need to construct the sequences of all possible intermediate states.This can be represented by vacuum Feynman diagrams.Each Feynman diagram is topologically distinct.For each Feynman diagram we move around the vertices, so as to consider all possible orderings and correspondingly all possible intermediate states.Without loss of generality, we can pick one of the vertices to be first, while considering different orderings of the other vertices.If we cut the Feynman diagram at this first

Figure 11 :
Figure 11: There is one tetrahedron diagram, which is related by symmetry to the tetrahedron diagram evaluated in Fig. 3.
(a) is the identity, Fig. 12 (b) is a rotation about vertex a, Fig. 12 (c) is a further rotation.Fig. 12 (d) is a mirror reflection of Fig. 12 (a) along the plane holding the 1 line fixed and bisecting the face bordered by lines 2, 3, 6.Fig. 12 (e) is a mirror reflection of Fig. 12 (b) holding the 2 line fixed, and Fig. 12 (f) is a mirror reflection of Fig. 12 (c) holding the 3 line fixed.Fig. 12(a) has vertices labeled a, b, c, d and edges labelled 1, 2,

Figure 12 :
Figure 12: The vertex orderings of the tetrahedra are: a) a, b, c, d; b) a, c, d, b; c) a, d, b, c; d) a, d, c, b; e) a, b, d, c; f) a, c, b, d .

Figure 13 : 2 + 2 +
Figure 13: (a) A one-loop diagram contributing to ⟨a p ⟩ (b) Higher order corrections to ⟨a p ⟩.The gray blob represents corrections to the propagator of ⟨a * p (t)a p (0)⟩.
(1,2,3,4,5,6ish, we can do a change of variables(1,2,3,4,5,6) → (4, 5, 1, 6, 2, 3) on the T b term, so as to rotate the b vertex into the a vertex to the greatest extent possible (meaning some of the arrows will be reversed; but we can't flip arrows under a change of variables.This gives back (C.1).Note that from this perspective, the second term in (C.1) makes sense, because if we had done the transformation taking us from a, b, c, d to b, d, a, c then the edges would have gone from 1, 2, 3, 4, 5, 6 to 3, 5, 6, 1, −2, 4 (see the table in Appendix.B), which is almost our transformation in the reverse direction, but with a minus sign for 5.