Classical Sub-subleading Soft Photon and Soft Graviton Theorems in Four Spacetime Dimensions

Classical soft photon and soft graviton theorems determine long wavelength electromagnetic and gravitational waveforms for a general classical scattering process in terms of the electric charges and asymptotic momenta of the ingoing and outgoing macroscopic objects. Performing Fourier transformation of the electromagnetic and gravitational waveforms in the frequency variable one finds electromagnetic and gravitational waveforms at late and early retarded time. Here extending the formalism developed in \cite{1912.06413}, we derive sub-subleading electromagnetic and gravitational waveforms which behave like $u^{-2}(\ln u)$ at early and late retarded time $u$ in four spacetime dimensions. We also have derived the sub-subleading soft photon theorem analyzing two loop amplitudes in scalar QED. Finally, we give the structure of leading non-analytic contribution to (sub)$^{n}$-leading classical soft photon and graviton theorems which behave like $u^{-n}(\ln u)^{n-1}$ for early and late retarded time $u$.

Since S-matrix for massless theory is IR divergent in four spacetime dimensions, analysis of soft theorem for loop amplitudes was known to be ambiguous [46][47][48][49]. But soft theorem relates two Smatrices, so one does not need to make individual S-matrices IR-finite, instead one can factor out the same IR divergent piece from both the S-matrices (if possible) and cancel the IR divergent piece from both sides in the soft theorem relation 1 . Using this prescription, soft photon and soft graviton theorems have been derived in [17] up to subleading order in soft momentum expansion. The soft factor at subleading order becomes logarithmic in the energy of external soft photon/graviton. This logarithmic soft factor turns out to be universal and one loop exact. In [17] the authors made an observation that in the Feynman diagrammatics of loop amplitude if one replaces the Feynman propagator for virtual photon/graviton by it's corresponding retarded propagator one gets only the classical soft factor at subleading order which is proportional to the classical electromagnetic/gravitational waveform 2 .
A direct evaluation of long wavelength electromagnetic and gravitational waveform has been performed in [1] for a general classical scattering process in the presence of long range electromagnetic and gravitational force. The strategy for deriving gravitational waveform followed in [1], is iteratively solving Einstein equation to find the metric fluctuation and geodesic equation to find the correction of asymptotic trajectories. Then performing Fourier transformation in the frequency of long wavelength gravitational waveform one finds the DC gravitational memory at leading order and u −1 gravitational tail memory at subleading order for retarded time u. Here in this work, we have given a systematic extension of the formalism developed in [1] for deriving classical soft photon and soft graviton theorems up to arbitrary order in soft momentum expansion. Using this prescription we have explicitly derived classical soft photon and soft graviton theorem at sub-subleading order which in turn gives us u −2 (ln u) tail memory after Fourier transformation. We have also shown how the same result could also be derived from two loop amplitudes following [17]. Finally we have given the structure of the leading non-analytic contribution of (sub) n -leading classical soft photon and soft graviton theorems in terms of some undetermined functions which reproduce u −n (ln u) n−1 tail memory.
We shall begin by stating our result of classical soft photon theorem, ignoring gravitational interaction.
Consider a scattering process in which M number of objects come in with momenta p 1 , p 2 , · · · , p M and charges q 1 , q 2 , · · · , q M undergo complicated interactions 3 in a large but finite region of spacetime and finally disperse to N number of outgoing objects with momenta p 1 , p 2 , · · · , p N and charges q 1 , q 2 , · · · , q N .
Let R be the region of spacetime of linear size L inside which all the complicated interactions happen and outside the region R only long range electromagnetic force acts between charged scattered objects. Now defining the scattering center to be some spacetime point well inside the region R and keeping the electromagnetic wave detector at distance R from the scattering center, we want to derive the 1 R component of the electromagnetic waveform at large retarded time |u| = |t − R| >> L. We are working in a unit where the speed of light c = 1. In the large u expansion the leading and subleading order electromagnetic waveforms are given in [1,17,21,[42][43][44][45] in four spacetime dimensions. At leading order the change of the electromagnetic waveform between u = −∞ and u = +∞ turns out to be constant (velocity kick memory for test charge) and at subleading order the early and late time electromagnetic waveforms go like u −1 as |u| → ∞ (acceleration memory for test charge). Here extending the analysis of [1], we derive the leading non-analytic contribution of sub-subleading electromagnetic waveforms at early and late time. The results are as follows: where, x ≡ Rn , n µ ≡ (1,n). (1.3) Above the expression for B (2) αµ q a , p a ; {q b }, {p b } takes form, with the expression for C α q a , p a ; {q b }, {p b } 4 : The expression for B (2) αµ q a , p a ; {q b }, {p b } has the same functional form with the arguments having charges and momenta for incoming objects and the particle sums run over M values. The above result is universal i.e. theory independent and determined in terms of four momenta and charges of scattered objects. On the other hand the O u −2 correction of the sub-subleading electromagnetic waveform turns out to be theory dependent. It also depends on the details of the scattering event and on the structures (via electromagnetic multipole moments) and spins of the objects participating in the scattering process.
The contribution of first two lines in the sub-subleading electromagnetic waveform of eq.(1.1) and (1.2) describes the effect of electromagnetic radiation due to late and early time acceleration of charged objects in the background leading electromagnetic field produced by a pair of two charged objects due to their asymptotic straight line trajectories. On the other hand the contribution B (2) αµ q a , p a ; {q b }, {p b } or B (2) αµ q a , p a ; {q b }, {p b } describes the effect of electromagnetic radiation due to the subleading correction to the straight line trajectories of charged objects. For massless charged particle scattering process, sub-subleading electromagnetic waveforms in eq.(1.1) and (1.2) vanish. Now consider the same scattering setup but take into account the effect of gravitational interaction.
For simplicity we shall consider the scattering objects are charge neutral so that outside the region R only long range gravitational force acts. So we are interested in deriving the 1 R component of the gravitational waveform at large retarded time |u| >> L. Let us define the deviation of metric from Minkowski metric as, In the large u expansion the leading and subleading order gravitational waveforms are given in [1, 21, 23-30, 33, 34, 51-55] in four spacetime dimensions. At leading order the change of the gravitational waveform 4 In the expression of Cα we can remove the terms proportional the paα as those terms ultimately cancel in the expression of B (2) αµ , but we kept those explicitly to make direct connection with Feynman diagram analysis as discussed in §4. between u = −∞ and u = +∞ turns out to be constant and at subleading order the early and late time gravitational waveforms go like u −1 as |u| → ∞. At sub-subleading order the leading non-analytic contribution of gravitational waveform has been conjectured in [1] which goes like u −2 (ln u) as |u| → ∞ at late and early time. Here extending the analysis of [1], we have given a proof of the conjectured subsubleading gravitational waveforms at early and late retarded time u = t − R + 2G ln R N b=1 p b · n, which has the following form, ∆ (2) e µν (t, R,n) and ∆ (2) e µν (t, R,n) The gravitational waveform given above at order O u −2 ln |u| is universal i.e. theory independent and determined in terms of the four momenta of scattered objects. On the other hand the un-determined O(u −2 ) gravitational memory depends on the spin angular momenta of the scattered objects as well as on the details of the scattering region R. As an observational prediction of our result, consider gravitational radiation from a binary black hole merger, which we can think of as a decay process with bound state of two black holes as the only initial state and one massive black hole and mass less gravitational radiation are the final states. For this process M = 1 and one final state object is massive, and rest are massless.
In [1], it is already discussed that with the order u −1 gravitational tail memory, the order u −2 ln |u| gravitational tail memory also vanishes for binary black hole merger problem. So the absence of both the tail memories for binary black hole merger is the prediction from Einstein gravity. But as discussed in [1] for supernova explosion, neutron star merger and hyper-velocity star production we can have nonvanishing finite gravitational tail memory.
The organization of the paper is as follows: In §2 after reviewing the derivation of leading and subleading classical soft photon theorem, we derive the sub-subleading classical soft photon theorem.
Then extending this formalism we give the structure of leading non-analytic contribution of (sub) n -leading classical soft photon theorem. In §3 after reviewing the derivation of leading and subleading classical soft graviton theorem, we derive the sub-subleading classical soft graviton theorem. In §4 we derive subsubleading soft photon theorem from the analysis of two loop amplitude and describe it's connection with sub-subleading classical soft photon theorem. In §5, first we give a Feynman diagrammatic understanding of the derivation of sub-subleading soft graviton theorem with the expected result. Then we give the structure of (sub) n -leading electromagnetic and gravitational waveform for a scattering process where both electromagnetic and gravitational long range forces act between the scattered objects. Finally we make some comments on the gravitational tail memory for spinning object scattering.
When this paper was largely complete, a paper [56] appeared where a derivation of sub-subleading electromagnetic waveform has been attempted by directly using position space Green's function. But the results of eq.(1.1) and eq.(1.2) has not been reproduced completely in eq.(61) of [56] after substituting the expression of f µ i from eq.(56) of the same paper. A direct comparison of f µ i expression in [56] with eq.(1.5) shows that the following contribution within the square bracket in the second line of RHS of eq.(1.5) has been missed in [56]  electromagnetic wave at distance R from the scattering center. Then in this setup we want to determine the 1 R component of the electromagnetic waveform for retarded time u >> L. This also translates to determining the radiative mode of gauge field with frequency ω in the range R −1 << ω << L −1 . As discussed in [1,[20][21][22]55] the Fourier transformation in time variable of the radiative mode of the gauge field is related to the Fourier transform of current density determined in terms of the asymptotic trajectory of scattered objects in the following way, where, and x = Rn withn being the unit vector along the direction of detector from the scattering center.
Since L must be bigger than the size of the objects involves in the scattering process and we are interested in determining electromagnetic waveform with wavelength larger than L, we can treat the objects involved in the scattering process as particles. This statement is true for determining the leading non-analytic contribution of (sub) n -leading electromagnetic waveform which turns out to be of order O ω n−1 (ln ω) n in small ω expansion. But if we want to determine the order O ω n−1 (ln ω) k contribution of (sub) n -leading electromagnetic waveform for k ≤ n − 1, the point particle assumption breaks down.

General setup
For performing the analysis for both ingoing and outgoing objects at one go we use the following compact notations: for a = 1, 2, · · · , M and 0 ≤ σ < ∞. Consider the asymptotic trajectory of a'th particle outside the region R, where {r a } denotes the set of points where particles' trajectory intersect the boundary of region R at σ = 0, as indicated in Fig.1. Now to find the asymptotic trajectory of the particles and the electromagnetic gauge field we have to solve the following two equations, where in the last equation above we have used Lorenz gauge condition ∂ µ A µ (x) = 0. The current density in terms of the asymptotic trajectories of the particles is given by, Here we are using the unit where speed of light c = 1 and the charges {q a } are in unit of √ α, where α is the fine structure constant. In this unit q 2 ω M is a dimensionless parameter with q being the parameter representing charges of scattered particles and M being the parameter representing masses of scattered particles 6 . So expansion of any quantity in power of ω is same as expansion of this quantity in powers of q 2 . Hence let us expand correction of the straight line trajectory in powers of q 2 in the following way, . Now with this correction of trajectory the Fourier transform of the current density can be expanded as, where ∆ (r) J µ (k) is of order O(q 2r+1 ). Similarly the gauge field also has an expansion, for r = 0, 1, 2, · · · . Here G r ( ) represents the momentum space retarded Greens function, Now from eq.(2.9), the expression of the Fourier transform of current density turns out to be, Here we want to extract the leading non-analytic contribution of ∆ (r) J µ (k) in the small ω expansion, which turns out to be of order O ω r−1 (lnω) r . For the analysis of leading and subleading order we will closely follow [1], but here we will be brief and organize the results in a different way which will be useful for the generalization to (sub) n -leading order.

Derivation of leading order electromagnetic waveform
The leading order current density follows form eq.(2.15), where for the convergence of the integral at σ = ∞ we need to replace ω → (ω + iη a ) where η a = ±1 for particle-a being outgoing/ingoing. This implies the following replacement in the exponential To extract the leading order contribution in ω expansion, we have to expand e −ik.ra = 1 − ik.r a + · · · and keep only the leading contribution. Now using the result of eq.(2.1) at leading order and performing Fourier transformation we recover the known leading electromagnetic memory, q a p aµ p a .n (2.17)

Derivation of subleading order electromagnetic waveform
From the leading order current density in eq.(2.16), the leading gauge field expression becomes, The expression for the leading order field strength produced due to the asymptotic straight line motion of particle-b is, For the leading order gauge field in eq.(2.18), leading correction to the straight line trajectory of particle-a satisfy the following equation of motion, After using the boundary conditions following from eq.(2.6) we get, From eq.(2.15), the expression of subleading order current density becomes, Now we will analyze each term within the square bracket above separately to extract contribution of order O(ln ω) following [1]. The first term of subleading current density in eq.(2.25) is, The self force on the trajectory of charged particles can be neglected up to the order we are working in [57]. where, Here in eq.(2.26), to get the second line of RHS from the first line we have first written Here also to get the last line of RHS from the second last expression, we approximated the integrand in the momentum region L −1 >> | µ | >> ω and the corresponding result is written using sign ignoring the O(ω 0 ) contribution. Now summing contributions from both the terms, total subleading current density 8 Under sign we are also neglecting terms of order O(ω k ln ω) for k = 1, 2, · · · which we can get by expanding exp(−ik.ra) = 1 − ik.ra + · · · from the second last line of eq.(2.26). But these terms will not affect the leading non-analytic contribution of ∆ (r) Jµ(k) for r = 2, 3, · · · . at order O(ln ω) becomes, where k µ = ωn µ with n = (1,n). Now using the relation (2.1) at subleading order and performing Fourier transformation in ω variable we get the following early and late time electromagnetic waveforms [1] after using the results of the integrations from eq.(C.10) and eq.(C.11) for n = 1,

Derivation of sub-subleading order electromagnetic waveform
Sub-subleading order current density from eq.(2.15) takes the following form, Now we will analyze each terms within the square bracket above separately to extract contribution of order O ω(ln ω) 2 , which turns out to be the leading non-analytic contribution at this order in the small ω expansion. The first term of sub-subleading current density in eq.(2.32) contributes to, Now after using the following integration result in the above expression Now the above expression contributes to order O ω(ln ω) 2 in the integration region Now in the above expression the d 4 2 integral can be evaluated first by performing 0 2 integration using contour prescription and then by performing the angular integral and finally doing the In the result of the d 4 2 integration, the 1 dependence turns out ln(| 1 |ω −1 ). Then following the same steps for d 4 1 , the omega dependence appears as ω × Clearly the origin of the half factor is due to the upper limit of d 4 2 integration being | 1 |, which comes from the momenta ordering | µ 1 | >> | µ 2 |. This suggests that we can also express the contribution of ∆ (2) J µ (k), approximated in the integration region L −1 >> | µ 1 | >> | µ 2 | >> ω, by both d 4 1 and d 4 2 integrals having lower limit ω and upper limit L −1 with including an overall multiplicative factor of 1/2. Following this equivalence prescription of organizing the two momentum integrations we get, Similarly the second term in the sub-subleading current density expression (2.32) contributes to, So using the same trick as described below eq.(2.36) we get, The third and fourth terms within the square bracket in eq.(2.32) also contribute to order O(ω(ln ω) 2 ) as explicitly evaluated in appendix-A. Summing the contributions of eq.(A.7) and eq.(A.9) we get, ∆ (2) J µ (k) + ∆ (2) J µ (k) with the expression for C α q a , p a ; {q b }, {p b } : Now summing over the contributions of eq.(2.37), (2.39) and (2.40), we get the following expression for the total sub-subleading current density at order O(ω(ln ω) 2 ), Now using the relation (2.1), the radiative component of the sub-subleading order electromagnetic waveform becomes, Above the expression of B αµ q a , p a ; {q b }, {p b } is same as the one given in eq.(2.41) but with the particle sums running over only outgoing objects and the expression of B (2) αµ q a , p a ; {q b }, {p b } takes the same functional form as in eq.(2.41) with the sums running over only incoming objects. Now taking Fourier transformation in ω variable and using the results of integration (C.10) and (C.11) for n = 2, we find the following late and early time sub-subleading electromagnetic waveform: with x = Rn, n = (1,n) and u = t − R.

Structure of (sub) n -leading electromagnetic waveform
In this subsection we derive the leading non-analytic contribution of (sub) n -leading electromagnetic waveform up to some un-determined gauge invariant function. Before analyzing the contribution of (sub) nleading current density, let us list some useful observations from the last three subsections. 2. For the analysis of the terms containing ∆ (2) Y a (σ) in the integrand of sub-subleading current density in eq.(2.32), we need to use the expression of subleading gauge field ∆ (1) A µ (x) as explicitly analyzed in appendix-A. We also found that sum of the contributions from those terms is gauge invariant by itself. So we expect that in the analysis of (sub) n -leading current density the part of the integrand containing ∆ (r) Y a (σ) for r ≥ 2 similarly contributes to order O ω n−1 (ln ω) n and the contribution will be gauge invariant by itself. But we will not be able to analyze those terms since for that we need the expressions of the gauge fields ∆ (r) A µ (x) for r = 1, ..., n − 1, which we don't have yet beyond r = 2. So we write the contribution from the analysis of those terms in terms of some undetermined function of momenta and charges, which can be derived by the method of induction.
The (sub) n -leading current density from eq.(2.15) turns out to be, perform integration by parts to move the σ derivative to the rest of the integrand after using the boundary conditions of eq.(2.22) we get, The first term within the square bracket contributes to, Now the above expression contributes to order O ω n−1 (ln ω) n in (n − 1)! numbers of integration region: where P denotes all possible permutations of i 's within the parentheses. In any of this (n − 1)! integration regions the leading contribution in ω expansion from the following integration becomes, Now substituting the result of the above integral in eq.(2.49), we find same contribution from all the integration regions. So the total contribution will be (n − 1)! times the contribution from one of the region, say L −1 >> 1 >> 2 >> · · · >> n−1 >> n >> ω. Now using the procedure described below eq.(2.36), the expression can also be written with all the n number of integrating momenta run from ω to L −1 by the cost of an overall 1 n! . Performing all these steps we find, Now if we compare the above expression with the second term in the sub-leading current density in eq.(2.28), we find that in the above expression the integrand contains an extra multiplicative factor of relative to the integrand of eq.(2.28). This suggest to include an extra factor of (1) J µ (k) following our first observation in the beginning of this subsection. So the above derivation gives a direct generalization of the first observation made in the beginning of this subsection. An analogous analysis for the second term within the square bracket of eq.(2.48) gives the following contribution at order O ω n−1 (ln ω) n , The contribution of ∆ rest (n) J µ (k) has not been evaluated explicitly due to the reason explained in the second observation in the beginning of this subsection but the general form follows from gauge invariance [58][59][60], with the property that B (n) } is anti-symmetric in µ and any one of the α r indices for r = 1, · · · , n − 1. We expect that B (n) α 1 α 2 ···α n−1 µ q a , p a ; {q b }, {p b } will be independent of the details of the scattering event as well as will be independent of the nature of short range forces acting between the scattered objects inside region R. After summing all the contributions, the O ω n−1 (ln ω) n contribution of ∆ (n) J µ (k) becomes, Now using the relation of eq.(2.1) at (sub) n -leading order, performing Fourier transformation in ω variable and using the results of Fourier transformations (C.10),(C.11) we get the following late and early time electromagnetic waveform, Above form the reality of electromagnetic waveform, we expect that (i) n B (n) is real and from momentum and charge power counting

Proof of classical soft graviton theorem
Consider the same scattering process as described in §2 but for simplicity consider the scattered objects are charge neutral. Now the region R is chosen such that all complicated interactions take place inside spacetime region R and outside only long range gravitational interaction is present as described in Fig.1.
We define the scattering center be a spacetime point well inside the region R of linear size L and place the detector of gravitational wave at distance R from the scattering center. Then our goal is to determine the 1 R component of gravitational waveform for retarded time u >> L. This translates to determining the radiative mode gravitational waveform with frequency ω in the range R −1 << ω << L −1 . Let us define the deviation of metric from Minkowski metric as, Now following [1,[20][21][22]55], the radiative component of the trace reversed metric fluctuation is related to the Fourier transform of the total energy-momentum tensor T µν (x) in the following way, where, and x = Rn , k µ = ωn µ = ω(1,n) withn being the unit vector along the direction of the detector from the scattering center.

General setup
We use the following compact notations for analyzing ingoing and outgoing particles' trajectories as done in §2.1, for a = 1, 2, · · · , M and 0 ≤ σ < ∞. Consider the asymptotic trajectory of a'th object outside the region R, where r µ a denotes the point on the boundary of the region R where particle-a enters or exits region R at σ = 0. Y µ a (σ) satisfies the following boundary conditions, Now we have to solve Einstein equation and geodesic equation 9 iteratively to determine gravitational waveform in terms of the scattering data.
where T Xµν is the matter energy-momentum tensor given by 10 , In de Donder gauge η ρµ ∂ ρ e µν = 0, the Einstein eq.(3.8) takes the following form, where T hµν (x) is the gravitational energy-momentum tensor defined as, Here we are following the unit convention c = 1 and in this unit GM ω is a dimensionless quantity, where ω is frequency of soft gravitational radiation and M is the parameter representing masses of scattered objects 11 . So expansion in power of ω is same as expansion in power of gravitational constant G. Let us expand the Fourier transform of total energy momentum tensor in a power series expansion of G, where ∆ (r) T µν (k) is of order (G) r in gravitational constant expansion. Similarly the correction to the straight line trajectory has an expansion in power of G, where ∆ (r) Y µ a (σ) is of order (G) r in gravitational constant expansion. This implies that the trace reversed metric fluctuation has the following expansion in power of G if we use (3.11), This implies ∆ (r) e µν (x) is of order (G) r+1 in gravitational constant expansion. Analogous to the analysis of electromagnetic waveform case here we expect that the leading non-analytic contribution of ∆ (n) T µν (k) should be of order ω n−1 (ln ω) n in small ω expansion and for extracting gravitational waveform at this order we can treat the scattered objects as non-spinning point particles 12 .
Fourier transform of matter energy momentum tensor in terms of corrected trajectory takes the following form, The gravitational energy-momentum tensor at order (G) r have the following kind of dependence on where ∂∂ above is just showing that each term in gravitational energy-momentum tensor has two derivatives on metric fluctuation 13 and the various {∆ (s) e} dependence is fixed from the requirement of ∆ (r) T hµν (x) should be of order (G) r . For r = 1 and r = 2 the expressions for energy-momentum tensor are explicitly given in appendix-D.
Though our final goal is to determine the sub-subleading gravitational waveform at order O(ω(ln ω) 2 ), to set up the stage we will go through the derivation of subleading order gravitational waveform briefly following [1] as some of the subleading order results are needed for the analysis of sub-subleading order energy-momentum tensor.

Derivation of leading order gravitational waveform
Fourier transform of the leading order matter energy-momentum tensor from eq.(3.17) takes form, Similar to §2.2, from the convergence of the integral at σ = ∞ for both ingoing and outgoing particles we need to replace ω → (ω + iη a ) where η a = ±1 for particle-a being outgoing/ingoing. This implies the following replacement in the exponential k · v a → (k · v a − i ) which we will use from now on. The gravitational energy momentum tensor at this order vanishes. Now using the relation (3.2) for the leading energy-momentum tensor and performing Fourier transformation in ω variable we get the expected DC gravitational memory [1,[23][24][25][26][27][28][29][30][31][32][33], Now from the leading order energy momentum tensor in eq. (3.19), the leading order metric fluctuation takes the form, The leading order Christoffel connection produced due to the asymptotic straight line trajectory of particle b is given by,

Derivation of subleading order gravitational waveform
For the leading order metric fluctuation in eq. (3.22), leading correction to the straight line trajectory of particle-a satisfies the following geodesic equation as follows from eq.(3.9), Now after using the boundary conditions (3.26) in eq.(3.24) and performing integration in σ variable we get, The expression for the subleading order matter energy-momentum tensor from eq.(3.17) takes form, For the analysis of this matter energy-momentum tensor we will be very brief and follow ref.
[1] mostly but analyze the terms within the square bracket separately and organize the results in a different way which will help us to get the structure to higher order matter energy-momentum tensor contribution in small ω expansion. In the first term within the square bracket of eq.
va−i )σ and then perform integration by parts to move the σ derivative to the rest of the integrand after using the boundary conditions from eq.(3.26), we get the following contribution, where in the last line above we have approximated the integrand in the integration range L −1 >> | µ | >> ω and written the approximated result using sign 15 . The expression for K cl gr is given below and the integration is explicitly evaluated in [1], Performing an analogous analysis for the second term within the square bracket of eq.(3.29) we get, where in the last line above we have approximated the integrand in the integration range L −1 >> | µ | >> ω under the sign . Similarly in this integration range the third term within the square bracket of eq.(3.29) contributes, Hence the order O(ln ω) contribution from subleading order matter energy-momentum tensor becomes, Fourier transform of subleading order gravitational energy-momentum tensor follows from (D.13), where, F µν,αβ,ρσ (k, ) Now following [1], ∆ (1) T hµν (k) contributes to ln ω in two integration regions. In the integration region R −1 << | µ | << ω with R being the distance of the detector from the scattering center, we can approximate F µν,αβ,ρσ (k, ) up to gauge equivalence as, Then the contribution from this integration region turns out to be, In the integration region ω << | µ | << L −1 , F µν,αβ,ρσ (k, ) can be approximated with the terms quadratic in . While evaluating this integral one will encounter a linearly divergent term in this region which can be identified with the leading order energy-momentum tensor for soft radiation from finite energy gravitational radiation. Since this kind of contribution is already taken care of inside the finite energy particle's sum we need to subtract this contribution by hand. The expression for this soft radiation energy-momentum tensor is Finally subtracting the contribution ∆ (0) T Rµν (k) from ∆ (1) T hµν (k), which is the contribution from ∆ (1) T hµν (k) in the integration region ω << | µ | << L −1 , we get, Now summing over eq.(3.34), eq.(3.40) and eq.(3.38) we get the total subleading energy-momentum tensor at order O(ln ω), Now to determine the radiative mode of subleading order gravitational waveform at late and early time we use the relation (3.2) at subleading order and then perform Fourier transformation in ω variable. Finally using the results of integrations from eq.(C.10) and (C.11) for n = 1 we get [1],

Derivation of sub-subleading order gravitational waveform
Here in this section we derive the order O(ω(ln ω) 2 ) coefficient of the sub-subleading gravitational waveform. This turns out to be the leading non-analytic contribution at this order in small ω expansion.

Analysis of sub-subleading matter energy-momentum tensor
The expression for sub-subleading order energy-momentum tensor follows from eq.(3.17), Contribution from the first term within the square bracket of eq.(3.44) takes form, Here to get second line of RHS from the first line we have first used | µ 1 | >> | µ 2 | >> ω. Now following the procedure described below eq.(2.36) we can also get the same contribution approximating the integrand in the integration region L −1 >> | µ 1 | >> ω , L −1 >> | µ 2 | >> ω but with an overall multiplicative factor of 1 2 coming from the momentum ordering. So the O ω(ln ω) 2 contribution becomes, Similarly the contribution from the second term within the square bracket of eq.(3.44) turns out to be, ∆ (2) T Xµν (k), the above expression also contributes to order O ω(ln ω) 2 in the integration region L −1 >> | µ 1 | >> | µ 2 | >> ω. So following the procedure described below eq.(2.36) we can analyze the two integrals in the integration region L −1 >> | µ 1 | >> ω , L −1 >> | µ 2 | >> ω but with an overall multiplicative factor of 1 2 coming from the momentum ordering. So the order O ω(ln ω) 2 contribution becomes, The contribution from the third term within the square bracket of eq.(3.44) can be easily read out from eq.(3.48) after interchanging µ ↔ ν and the order O ω(ln ω) 2 contribution is, Contribution from the fourth term within the square bracket of eq.(3.17) takes the following form after performing the σ integrations, Analyzing all possible integration regions we find that the above expression does not contribute at order O(ω(ln ω) 2 ). The expression ∆ (2) T Xµν (k) starts contributing at order O(ω ln ω) in the integration region It also contributes to order O(ω 0 ) but since order O(ω 0 ) is already undetermined from the analysis of subleading order and does not contribute to memory, we are ingoing that contribution 16 .
Analysis of fifth, sixth and seventh terms within the square bracket of eq.(3.44) is little involved as we need to find out subleading correction ∆ (2) Y µ a (σ) to the straight line trajectory after knowing subleading order gravitational waveform ∆ (1) e µν (x). A complete analysis has been performed in appendix-B and those terms also contribute to order O(ω(ln ω) 2 ).

Analysis of sub-subleading gravitational energy-momentum tensor
Now we want to analyze the sub-subleading order gravitational energy-momentum tensor to extract O(ω(ln ω) 2 ) contribution. From eq.(3.18) if we substitute r = 2, we find that there are two kinds of contributions to sub-subleading order gravitational energy-momentum tensor. One is cubic in ∆ (0) h µν which we denote by ∆ (2) T hµν 1 (x) and the other contribution is quadratic in metric fluctuation with one of them is ∆ (0) h µν and the other is ∆ (1) h µν which we denote by ∆ (2) T hµν 2 (x). The Fourier transform of the first part of the sub-subleading gravitational energy-momentum tensor follows from eq.(D.15), 16 One point to note that order O(ω 0 ) contribution always dominates over order O(ω(ln ω) 2 ) contribution. But if we are only interested to extract O(ω(ln ω) 2 ) coefficient unambiguously, we can do so by operating ∂ω on the integral expression and extract the coefficient of (ln ω) 2 .
where G µν αβ,ρσ,γδ ( 1 , 2 , k) can be read off from eq.(D.14) and the expression is quadratic in momenta. Since T (3)hµν (x) of eq.(D.14) contains two derivatives on one or two different {h τ κ }, we can always choose h γδ to be the metric fluctuation where no derivative operates, as explained in appendix-D. With this choice the momentum dependence of G µν αβ,ρσ,γδ ( 1 , 2 , k) turns out to be, A detail power counting analysis shows that ∆ (2) T hµν The Fourier transformation of ∆ (2) T hµν 2 (x) takes the following form, with the expression for ∆ (1) T µν ( 1 ) follows from eq.(3.29) and eq.(3.35) : where the expression for F µν,αβ,ρσ is given in eq.(3.36) and the expression for for E µν is the following, Here our goal will be to analyze the expression (3.54) in all possible integration regions after substituting the expressions of {∆ (1) h αβ } and {∆ (0) h αβ } to extract the order O ω(ln ω) 2 contribution. We only found the following three integration regions where eq.(3.54) contributes to order O ω(ln ω) 2 : In the integration region L −1 >> | µ 2 | >> ω >> | µ 1 | >> R −1 , approximating the expression in eq.(3.54) and performing the integration over 2 analogous to §3.3 from the first term within the square bracket we get, Now using the result of the integration: the contribution of U 1 turns out to be, In the region of integration integration ω >> | µ 2 | >> | µ 1 | >> R −1 , approximating the expression in eq.(3.54) and performing the integration over 2 analogous to §3.3 we find, Now using the result of the following integral: we get, For analyzing the expression (3.54) in the integration region L −1 >> | µ 2 | >> | µ 1 | >> ω, we need to substitute First let us analyze the contribution of eq.(3.54) with G * r ( 1 )G r (k − 1 ) = G r (− 1 )G r (k − 1 ). It turns out that the following part of ∆ (2) T hµν 2 (k) contributes to order O ω(ln ω) 2 in the integration region In the above expression ∆ ( 1 1 ) F µν αβ,ρσ (k, 1 ) represents the order O( 1 1 ) contribution of F µν αβ,ρσ (k, 1 ) and ∆ (k 1 ) F µν αβ,ρσ (k, 1 ) represents the order O(k 1 ) contribution of F µν αβ,ρσ (k, 1 ), which one can read out from the expression (3.36). Now it turns out that if we pick up the leading contributions from all the square brackets above in the integration region L −1 >> | µ 2 | >> | µ 1 | >> ω, the resulting expression does not contribute to order O ω(ln ω) 2 . On the other hand if we take subleading contribution from one of the square bracket above and leading contributions from the rest of the square brackets that can contribute to order O ω(ln ω) 2 . But for analyzing those contributions we have to use the following results which we get using integration by parts, Now after using this integration by parts results in U 3 and evaluating the integrals in the integration (2) T Xµν (k) − ∆ (2) T Xµν (k) We have to redefine the dummy indices for particle sums a, b, c and organize the huge expression coming from eq.(3.66) after performing the integrations, to bring the result in the compact form written above.
Now we also have to analyze the contribution of ∆ (2) T hµν as follows from eq.(3.65). With this substituation ∆ (2) T hµν 2 (k) can contribute to some linearly divergent term at order O(ω 0 L −1 ) in the integration region where integrating momenta are large compare to ω. We expect this contribution will represent the subleading soft gravitational radiation from hard gravitational radiation which needs to be subtracted by hand as it's effect is already included in the sum over hard particles a, b, c 17 . After this subtraction we find the following integral structure from ∆ (2) T hµν 2 (k), which can potencially contribute to O(ω(ln ω) 2 ), where f, g, h are functions of p a , p b , p c , n and the structure of the integral is written after performing integration over 2 . Now one can see that the part of the integrand inside the square bracket is symmetric  17 It would be interesting to independently compute subleading soft gravitational radiation from real hard gravitons at subleading order and compare with the O(ω 0 L −1 ) terms which we need to subtract here analogous to the leading order case as given in eq.(3.39).

Sub-subleading order gravitational waveform
Now using the relation in eq.(3.2), the sub-subleading order gravitational waveform in frequency space becomes, In the LHS of the above expression we have included a phase factor exp{−2iG N a=1 p a .k ln R} which corresponds to the time delay due to gravitational drag force on the emitted gravitational wave as derived in [17]. The expansion of this extra phase factor up to order ω 2 has already appeared in the expression (3.71). Now performing Fourier transformation in ω variable and using the integration results of eq.(C.11)(C.10)(C.16) we get the following late and early time sub-subleading gravitational waveforms, for u → +∞, (3.73) and p b · n is the retarded time. This proves the conjecture made in [1]. For a scattering event with massless external particles this result has been tested in [76,77].

Sub-subleading soft photon theorem from two loop amplitudes
Here in this section we briefly discuss how we can derive sub-subleading soft photon theorem by analyzing two loop amplitudes in a theory of scalar QED following the analysis of [17]. This will solve two purposes: Though in four spacetime dimensions S-matrix for this theory is IR divergent, still we can factor out the same IR-divergent piece from the S-matrix with external photon and S-matrix without external photon and can derive the soft photon factor relating the IR finite parts of the two S-matrices [17]. Let Γ (M +N,1) denote the full perturbative scattering amplitude of M + N number of external charged particles with momenta {p a }, charges {q a } and one external outgoing photon with momentum k. Similarly Γ (M +N ) denotes the full perturbative scattering amplitude with the same M + N number of external charged particles but without any external photon. Following [78], let us decompose the photon propagator attached with two complex scalar lines with momenta p a and p b in the following way, with the photon momentum runs from particle-b to particle-a. We are using the convention that all the external particles are outgoing and some of them can be made incoming by flipping the sign of their charges and momenta. The expression for K µν (ab) ( , p a , p b ) and G µν (ab) ( , p a , p b ) for a = b are, For a = b we do not need to make any decomposition. The propagator with K (ab) numerator is named as K-photon propagator and the propagator with G (ab) numerator is named as G-photon propagator. It turns out that if we compute the amplitudes Γ (M +N,1) and Γ (M +N ) after the K-G decomposition of the photon propagators described above we get [78],  In principle, in the LHS of the above relation we should also include the contribution of the two loop Feynman diagrams which contain one or both virtual photons whose both ends are connected to a single scalar line. But a detailed analysis shows that those diagrams don't contribute to order O(ω(ln ω) 2 ). Now one of the important observation made in [17] is that if we replace the G-photon Feynman propagators by G-photon retarded propagators then analyzing loop amplitudes with retarded G-photon  G,2−loop and use the algorithm of eq.(4.6) we expect the following relation will hold, So let us first re-derive our classical result of §2.4 analyzing Γ (M +N,1) G,2−loop , under the replacement rule described above.
• Contribution of diagram-(b) of Fig.2 reproduces the expression of ∆ (2) (2) J µ (k) in eq.(2.39) as well as the first three lines of the RHS of ∆ • On the other hand, sum over the diagrams of Fig.3 reproduces sum over the contributions ∆ (1) via the relation (4.7).
The one to one mapping between Feynman diagrams and classical computation described above is only true when the integrals are approximated in the specified regions of integrations from where the O(ω(ln ω) 2 ) contribution comes.
If we analyze the diagrams in Fig.2 and Fig.3 using G-photon Feynman propagators then using the algorithm in eq.(4.6) we find the full(classical +quantum) sub-subleading soft photon factor, which turn out to be: and where, . (4.12) The first and second lines of the RHS of sub-subleading soft photon factor expression (4.8), are gauge invariant separately. In principle the last line of eq.(4.8) can also be written in the following form: Hence, the origin of the above expression can be thought of as some effective coupling of electromagnetic field strength F µα (k) = i(k µ ε α − k α ε µ ) to the finite energy charged particles. Though the leading nonanalytic contribution to the sub-subleading soft photon factor given in eq.(4.8) is derived by analyzing two loop amplitudes in scalar QED, we expect this result is universal and 2-loop exact. Similarly if we generalize this analysis to n-loop amplitudes, we can in principle able to derive the leading non-analytic contribution of (sub) n -leading soft photon theorem whose classical limit will reproduce the result of eq.(2.54).

Outlook and speculations: further understanding of electromagnetic and gravitational waveforms
In this section the results we are giving are not rigorously derived. But under some specified assumptions we expect that the results are correct, but need future investigation.

Feynman diagrammatic understanding of the sub-subleading soft graviton theorem
Emboldened by our understanding of classical sub-subleading soft photon factor from the analysis of two loop amplitude in §4 with the replacement of Feynman propagators by retarded propagators, we expect that a similar analysis of two loop amplitudes for scalar coupled to gravity will reproduce classical subsubleading soft graviton factor. But as described in [17], for the derivation of soft graviton theorem we need the following four extra ingredients relative to the derivation of soft photon theorem.
1. We have to perform K-G decomposition of the graviton propagator, whose both ends are attached to same scalar line.
2. We do not need to perform any K-G decomposition for the graviton propagators whose one or both ends are connected to any graviton line. residual comes from the residual non-factorized contribution as discussed in the third point above.
• Contribution of Γ (M +N,1) 3,4−graviton will come from two kinds of diagrams. One kind of diagrams contain one or two three graviton vertices which can contribute to order O(ω(ln ω) 2 ) as we have seen from classical analysis. On the other hand from the classical analysis it is suggestive that the other kind of diagrams containing four graviton vertex do not contribute to order O(ω(ln ω) 2 ).
We expect that with this understanding one can extract sub-subleading soft graviton factor using the algorithm (5.1). Our expectation for the sub-subleading soft graviton factor is: where, 5.2 Structure of (sub) n -leading electromagnetic and gravitational waveforms in presence of both long range electromagnetic and gravitational forces In this section we take care of both electromagnetic and gravitational long range forces between the charged scattered objects outside the region R to give the structure of electromagnetic and gravitational waveforms at late and early retarded time. The trajectories outside the region R satisfies, Here we have to expand the asymptotic trajectories of scattered objects in the same way as in eq.(2.10), considering q 2 and G being the expansion parameters with equal weight.
Structure of (sub) n -leading electromagnetic waveform To find the electromagnetic waveform we have to solve the following form of Maxwell equation in the background metric produces by other particles, Using Lorentz gauge condition η µν ∂ µ A ν = 0 this can also be written as, where the gravitational contribution to current density is defined as, Now we have to expand the total current density in a power series of q 2 and G with both having the same importance as expansion parameter. Analogous to §2.5, here we expect the order O(ω n−1 (ln ω) n ) contribution to (sub) n -leading current density will take the following form 18 , ln(ω + i η a ) r n α 1 n α 2 · · · n α r−1 B (r),α 1 α 2 ···α r−1 µ q a , p a (5.9) where the undetermined function B (r),α 1 α 2 ···α r−1 µ q a , p a is antisymmetric under µ and any α exchange for = 1, 2, ..., r − 1. We expect that B (r) will be independent of the details of the scattering event inside 18 A brief sketch of the derivation of this result is given in appendix-E. the region R and will only depend on the scattering data. Now using the following relation (a modification of (2.1)) and performing Fourier transformation in ω variable after using the results of eq.(C.10),(C.11) and (C. 16) we find the following late and early time (sub) n -leading electromagnetic waveforms for n ≥ 2, n α 1 n α 2 · · · n α r−1 B (r),α 1 α 2 ···α r−1 µ − q a , −p a + O u −n (ln |u|) n−2 for u → +∞ (5.11) and where the retarded time u is given by It would be interesting to find out the expressions of {B (r) } generalising our analysis to arbitrarily higher order.
Structure of (sub) n -leading gravitational waveform Analysis of (sub) n -leading gravitational waveform for arbitrary value of n will be very complicated as higher order gravitational energy-momentum tensor will carry lots of terms. But we can still give some structure of gravitational waveform compiling the results of [58][59][60] and the observation made in [17].
In the classical analysis we expect the structure of order O(ω n−1 (ln ω) n ) contribution from the (sub) nleading matter and gravitational energy-momentum tensor for n ≥ 3 will take the following form 19 , × n α 1 n α 2 · · · n α r−1 C (r), µνα 1 ···α r−1 (q a , p a ) (5.13) where C (r), µνα 1 ···α r−1 (q a , p a ) is anti-symmetric under µ ↔ α i exchange as well as ν ↔ α j for i, j = 1, 2, · · · , r − 1. The functional behaviour of C (r) (q a , p a ) has to determine by explicit computations following eq.(E.10), but we expect these functions will only depend on the scattering data. Using the relation of eq.(3.2) and performing Fourier transformation in ω variable after using the results of Fourier transformations given in eq.(C.10),(C. 16) and (C.11) we get the following late and early time (sub) n -leading gravitational waveforms for n ≥ 3, q a q e p 2 e p 2 a p e .n − p a .p e p a .n 2 p a .p e p 2 a p 2 e p e .n − (p a .p e ) 3 p e .n − 1 2 p 2 a (p 2 e ) 2 p a .n + 1 4π q a q e p 2 e p 2 a p e .n − p a .p e p a .n and

Comments on gravitational tail memory for spinning object scattering
In the analysis of §3 we made a statement that internal structures or spins of scattered objects will not affect the leading non-analytic contribution to (sub) n -leading gravitational waveform. In this section we will try to present some arguments in support of this and also try to explain why determination of spin dependent gravitational memory is hard and ambiguous in our formalism.
We know that matter energy-momentum tensor has the following derivative expansion [22,68,79], where Σ a represents the spin of object-a and "· · · " contains terms having two or more derivatives operating on the delta function. The terms inside "· · · " carry the information about the internal structure of the scattered objects. For the spinning objects the geodesic equation is modified in the following way [80][81][82], where D dσ denotes covariant derivative along the world line. In the second equation above, we can substitute the leading large σ behaviour for Christoffel connection which goes like Γ µ νρ ∼ 1 σ 2 for large σ.
This in some sense gives conservation relation between different spin components, So, to analyze the leading order spin dependence of matter energy-momentum tensor we can treat the spins of scattered objects {Σ a } to be constant. Now at the leading order considering the trajectory of object-a being X µ a (σ) = r µ a + v µ a σ, the Fourier transformation of energy-momentum tensor becomes, Using the relation (3.2), we find that the leading spin dependent contribution to e µν (ω, R,n) appear at order O(ω 0 ) , which does not contribute to memory.
Solving linearized Einstein equation the leading order metric fluctuation takes form, with constraint on the spin p α a Σ a,αβ = 0. Now we can find out the subleading correction to the straight line trajectories solving geodesic equation of scattered objects in the background of the above metric analogous to §3.3. Here the subleading contribution to the Fourier transform of matter energy-momentum tensor becomes, Now if we try to analyze the first term within the square bracket above, we find the following spin dependent contribution, In the region of integration L −1 >> | µ | >> ω, the leading contribution turns out to be at order Similar analysis of the other terms in ∆ T Xµν (k) shows that those also contribute to O(ω ln ω). We can also analyze the subleading gravitational energy-momentum tensor, A detailed analysis shows that this will also contribute to O(ω 0 ) and O(ω ln ω) similar to ∆ (1) T Xµν (k).
So from this naive analysis we find that if the scattered objects carry spin then that will affect the order from the other non-universal contributions at the same order.
Though we have not been able to extract spin dependent gravitational memory unambiguously in this subsection but the analysis shows that the spins of the scattered objects do not affect our sub-subleading order gravitational waveform at order O(ω(ln ω) 2 ). Analogous analysis of the (sub) n -leading gravitational waveform shows that the spins of the scattered objects start contributing to gravitational tail memory at order O u −n (ln u) n−2 . Since the internal structure dependence comes with more powers of momenta in the Fourier transform of matter energy-momentum tensor follows from eq.(5.16), those start affecting the (sub) n -leading gravitational memory at order O u −n (ln u) m for m ≤ n − 3. A Analysis of the terms containing ∆ (2) Y a (σ) in the expression of ∆ (2) J(k) Before we proceed for analyzing the third and fourth terms within the square bracket of eq.(2.32), let us find out the expression of Subleading current density for b'th particle motion from eq.(2.25), With this form of subleading current density, the subleading field strength for b'th particle's motion takes the following form, The order O(q 4 ) correction to the straight line trajectory satisfies the following eq. of motion, Now using the boundary conditions follows from eq.(2.6): we get, The contribution to sub-subleading current due to third term in the square bracket of eq.(2.32), Analyzing the above expression in all the integration regions we find that only the last six lines above contribution from the last six lines becomes, Similarly the contribution to sub-subleading current due to fourth term in the square bracket of eq.(2.32) takes form, Analogously, analyzing the above expression in all the integration regions we find that only the last six lines contribute to order O(ω(ln ω) 2 ) in the integration region L −1 >> | µ 2 | >> | µ 1 | >> ω. Now using the prescription described below eq.(2.36) we find the following order O(ω(ln ω) 2 ) contribution, where the expression for F µν,αβ,ρσ is given in eq.(3.36) and the expression for for E µν is given as, Corresponding to this subleading energy-momentum tensor for b'th particle, the subleading metric fluctuation becomes, The order O(G 2 ) correction to the straight line trajectory of particle-a satisfies the following equation, Now after extracting the order O(G 2 ) contribution from the RHS above equation and using the boundary conditions follows from eq.(3.7): we get,

C Fourier transforms for deriving early and late time waveforms
Here we are interested to determine the Fourier transformation of following three functions,   where u ≡ t − φ. For u > 0 from the convergence of the Fourier integrals we have to close the contour in the lower half plane of complex ω. Similarly for u < 0 we have to to close the contour in the upper half plane of complex ω. Now for standard principle value of complex logarithms ln(ω ± i ), the branch cut singularities start at ω = ∓i and extend to ω = −∞ ∓ i . This implies F n,0 (t, x) = 0 for u < 0 and F 0,n (t, x) = 0 for u > 0. Now using the discontinuity of principal valued logarithm, defined as ln z = ln |z| + iArg(z) for −π < Arg(z) ≤ π, we get ln(ω + i ) = ln(ω − i ) + 2πiH(−ω), where H is the Heaviside step function. Substituting ln(ω + i ) = ln(ω − i ) + 2πiH(−ω) in eq.(C.4) we get, The most singular term in the integrand contributes,    Then from the last two expressions we find the following relation,

D Gravitational energy-momentum tensor
We consider deviation of the metric h µν from Minkowski metric given as g µν = η µν +2h µν and h µν satisfies harmonic gauge ∂ µ h µν − 1 2 ∂ ν h = 0. The definitions for Christoffel symbol, Riemann tensor and Ricci scalar are following: To compute various quantities below in terms of linear perturbation we follow the references [83][84][85].
E Brief sketch of the derivation of the results in §5.2 Here in this appendix we briefly discuss the origin of different terms in eq.(5.9) and eq.(5.13) for the scattering event described in §5.2 in presence of long range electromagnetic and gravitational interactions.
In presence of both long range electromagnetic and gravitational interactions, the Fourier transformation of (sub) n -leading order matter current density given in eq.(2.47) contributes to order O ω n−1 (ln ω) n in the integration regions where all the integrating momenta are larger than ω. Following the analysis of ln{(ω + i η a )L} r × n α 1 n α 2 · · · n α r−1 C (r), µνα 1 ···α r−1 (q a , p a ) (E.9) where C (r), µνα 1 ···α r−1 (q a , p a ) is anti-symmetric under µ ↔ α i as well as ν ↔ α j exchange for i, j = 1, 2, · · · , r − 1. The functional behaviour of C (r) has to determine by explicit computations, but we expect these functions will only depend on the scattering data. On the other hand in the integration region L −1 >> | µ 1 | >> ω, analysis of the (sub) n -leading gravitational energy-momentum tensor of eq.(E.8) is very much involved. But from the gauge invariance requirement of the total (sub) n -leading energy-momentum tensor, we expect the following structure of contribution at order O(ω n−1 (ln ω) n ) in the integration region where the integrating momenta are large compare to ω, ln{(ω + i η a )L} n n α 1 n α 2 · · · n α n−1 C (n), µνα 1 ···α n−1 (q a , p a ) (E.10) By analyzing the LHS of the above expression for n ≥ 3 in the specified integration region we can determine the unknown function C (n) (q a , p a ). Now summing over the contributions of eq.(E.6), eq.(E.9) and eq.(E.10) we get the expression of eq.(5.13).