Universality of loop corrected soft theorems in 4d

In [1], logarithmic correction to subleading soft photon and soft graviton theorems have been derived in four spacetime dimensions from the ratio of IR-finite S-matrices. This has been achieved after factoring out IR-divergent components from the traditional electromagnetic and gravitational S-matrices using Grammer-Yennie prescription. Although the loop corrected subleading soft theorems are derived from one-loop scattering amplitudes involving scalar particles in a minimally coupled theory with scalar contact interaction, it has been conjectured that the soft factors are universal (theory independent) and one-loop exact (don’t receive corrections from higher loops). This paper extends the analysis conducted in [1] to encompass general spinning particle scattering with non-minimal couplings permitted by gauge invariance and general coordinate invariance. By re-deriving the ln ω soft factors in this generic setup, we establish their universal nature. Furthermore, we summarize the results of loop corrected soft photon and graviton theorems up to sub-subleading order, which follows from the analysis of one and two loop QED and quantum gravity S-matrices. While the classical versions of these soft factors have already been derived in the literature, we put forth conjectures regarding the quantum soft factors and outline potential strategies for their derivation.


Introduction and Result
The soft theorem examines the infrared properties of a scattering amplitude involving a low (soft) momentum photon or graviton, in addition to other asymptotic particles. It establishes a relationship between this amplitude and the one without the low momentum photon or graviton. In a series of papers [2][3][4][5][6], it has been established that tree level soft photon and soft graviton theorems in four spacetime dimensions (D = 4) are just the manifestation of gauge invariance and general coordinate invariance at the scattering amplitude level. Soft factorisation alone does not provide profound insights into the ultraviolet completion of QED or quantum gravity theory, nor does it impose additional constraints on the quantum theory beyond what has already been achieved by gauge invariance and general coordinate invariance. 1 Instead, given an effective field theory (EFT) action with potential non-minimal interactions permitted by gauge invariance or general coordinate invariance, one can systematically compute the non-universal soft factors up to a certain order in the soft momentum expansion [3,5,8]. In the past, there were challenges in obtaining loop corrections to the subleading soft photon and graviton theorems in D = 4 due to the presence of infrared divergence in traditional scattering amplitudes [9,10]. However, this issue has been successfully resolved in [1] by directly working in D = 4 and carefully analyzing the possible non-analytic structures around ω = 0. In this work, the subleading soft photon and soft graviton theorems have been derived at the one-loop level, taking into account both electromagnetic and gravitational interactions. Interestingly, it has been observed that the subleading soft factors emerge at an order ln ω, where ω represents the energy of the soft photon or graviton. The loop corrected subleading soft factor is dominant compare the tree level subleading soft factor, which is of order ω 0 , as the energy approaches zero (ω → 0). The existence of the ln ω soft graviton theorem has been confirmed in [11,12] with perfect agreement with the result of [1] in the massless limit.
The soft graviton theorem results offer an intriguing application in deriving low-frequency gravitational wave forms and gravitational memory for astrophysical scattering events from their classical limit [13][14][15][16][17][18][19][20][21][22][23]. In a typical classical gravitational scattering scenario, one provides initial scattering data such as masses, velocities, sizes, intrinsic angular momenta, and impact parameters of the scattered objects, along with the specified interaction among them. The goal is to determine the gravitational waveform as an output. However, the classical limit of the universal soft graviton theorem directly provides the low-frequency gravitational waveform in terms of both the initial and final scattering data, regardless of the knowledge of the interaction involved in the scattering process. This suggests a novel approach for deriving low-frequency and late-time gravitational waveforms by directly studying classical gravitational scattering processes with both initial and final scattering data, known as the classical soft graviton theorem. This approach has been successfully pursued in [15,16,21,23]. The derivation of the classical soft graviton theorem readily extends to higher orders in the low-frequency expansion of the gravitational waveform, and numerous higher-order terms have been derived. However, deriving their quantum counterparts from the analysis of scattering amplitudes proves to be challenging in general, as discussed in the main body of the paper.
Consider a gravitational scattering amplitude involving N number of finite energy particles (hard particles) with momenta, spins and polarizations {p i , Σ i , ϵ i } for i = 1, 2, · · · , N 1 We would like to emphasize that the Ward identities relating these soft theorems to asymptotic symmetries also do not provide any additional constraints on the quantum theory of electromagnetism or gravity beyond what has already been achieved by gauge invariance and general coordinate invariance [7]. Instead, the Ward identities are simply the manifestation of the equations of motion for low-frequency photons or graviton fields at the level of scattering amplitudes. Where the soft charge is represented as an integral over the radiation mode of gauge or graviton fields and the hard charge is an integral over the inverse propagator operating on the current or stress tensor associated with finite energy scattered particles. The asymptotic symmetry parameters (large gauge transformation and asymptotic radial modes of bulk diffeomorphism) on the celestial sphere are just smearing functions that appears in both the integrands of the soft and hard charge expressions.
(N ≥ 4) and one low-energy (soft) graviton 2 with momentum and polarization k, ε, and denote this scattering amplitude by A (N +1) {ϵ i , p i , Σ i }, ε, k . We are following the convention that all the particles are ingoing, so if some of the particles are outgoing then we have to flip the sign of the four momenta for those particles. We will treat the soft graviton as outgoing with energy ω, so according to our convention k µ = −ωn µ with n µ being the null vector whose spatial part denotes the direction of soft graviton emission. Now the soft expansion of this (N + 1) particle amplitude takes the following form where A (N )β (i) (p i ) represents the i'th particle polarisation (ϵ iβ ) stripped N -particle amplitude A (N ) {ϵ i , p i , Σ i } , which is defined by the following relation (1.2) In (1.1), the expression of tree level "soft factor" 3 for single soft graviton emission reads [2,3,8,[24][25][26][27][28][29][30][31][32][33] (S gr In D = 4 analyzing tree level scattering amplitudes for effective field theory the nonuniversal sub-subleading soft graviton factor at order ω has also been derived in [3,8]. Reference [3] explicitly evaluated the non-universal contribution to the sub-subleading soft factor in terms of the non-minimal coupling of two finite energy fields to a soft graviton field through the Riemann tensor, and the general structure of the three-point 1PI (oneparticle irreducible) vertex involving two hard particles and a soft graviton. By extending the analysis of [3], it becomes evident that a complete soft factorization is not achievable at order ω n for n ≥ 2 in a generic theory of quantum gravity involving all possible higher derivative corrections allowed by general covariance. However, a partial soft factorization has been accomplished by enforcing linearized gauge invariance of the (N + 1)-particle amplitude in [6,34]. The generalization of the tree-level soft factor S gr tree for multiple soft graviton emissions up to subleading order can be found in [4,5,[35][36][37][38][39][40][41][42][43].
The one-loop contribution to the "soft factor" for single soft graviton emission in (1.1) reads and (1. 6) In the expressions (1.5) and (1.6), η j = +1 if j'th particle is ingoing and η j = −1 if j'th particle is outgoing. Under the sign ≃ we only keep the logarithmic contributions after performing the integrations following [1]. The upper limit of the integration Λ in (1.5) represents the order of the energy of hard particles and the lower limit of the integration R −1 in (1.6) represents the order of the energy resolution of the detector. The first line of (1.4) contains the O(ln ω) soft factor which has been derived in [1] as an one-loop exact result, analyzing one-loop gravitational S-matrices in the theory of scalar coupled to gravity.
In [1] a correction to O(ln ω) soft graviton factor due to electromagnetic interaction has also been derived when the scattering particles carry some electric charges as well. In this article we re-derive the O(ln ω) soft factor for single soft graviton emission in a generic theory of quantum gravity for scattering of particles with arbitrary mass and spin. This investigation will demonstrate the universal (independent of theory) nature of the O(ln ω) soft factor, while also extending the infrared divergence factorization prescription proposed in [1] to encompass a broad range of quantum gravity theories. In this paper we also conjecture the order ω ln ω soft factor given in the second and third lines of (1.4) which is derivable from the analysis of one-loop scattering amplitude for the scattering generic spinning particles in a generic theory of quantum gravity extending the analysis of this article. The classical limit of this O(ω ln ω) soft factor has already been derived in [23] in the name of spindependent classical soft graviton theorem which provides evidence on the correctness of the conjecture. Note that the tree level subleading soft theorem result at order ω 0 in (1.3) is not universal as it receives correction at one-loop order, which is expected to be dependent on the theory as well as the value of detector resolution (IR regulator) [9,10]. The two-loop contribution to the "soft factor" for single soft graviton emission in (1.1) reads (1.7) The above result was conjectured in [21] as a two-loop exact result and can be obtained by analyzing two-loop amplitudes using the same methodology being developed in this paper.
The classical limit of this order, denoted as O ω(ln ω) 2 , was derived under the name of the sub-subleading classical soft graviton theorem in [21]. This classical derivation offers substantial evidence supporting the validity of the above two-loop soft factor. Note that the O(ω ln ω) soft factor in the second and third lines of (1.4) at one-loop receives a correction at two-loop order. This correction is expected to depend on the theory of quantum gravity and the value of detector resolution (IR regulator), hence non-universal. From the analysis of the n-loop gravitational S-matrix, it is expected that the new leading non-analytic soft graviton factor, as the frequency ω approaches zero, behaves like ω n−1 (ln ω) n relating it to tree level N -particle amplitude. The general structure of the order ω n−1 (ln ω) n soft graviton theorem is provided in [21]. Note that the "tree", "1-loop" , ... subscripts in the soft factors in the expression (1.1) only specifies the first appearance of the soft factors in the perturbative analysis of the (N + 1)-particle amplitudes at that order (tree or n-loop) and it multiplies to the corresponding tree-level N -particle amplitudes. But they also appears in the analysis of higher loop amplitudes as well. For example S gr tree also appears as a soft factor in the analysis of (N + 1)-particle n-loop amplitude relating it to the N -particle n-loop amplitude for all n ≥ 1. Similarly S gr 1-loop also appears as a soft factor in the analysis of (N + 1)-particle n-loop amplitude relating it to the N -particle (n − 1)-loop amplitude for all n ≥ 2, and S gr 2-loop also appears as a soft factor in the analysis of (N + 1)-particle n-loop amplitude relating it to the N -particle (n − 2)-loop amplitude for all n ≥ 3. These observations also apply to the soft photon theorem results provided below. Now let us consider the same setup of scattering but turn off the gravitational interaction and turn on electromagnetic interaction between charged particles. We consider the finite energy scattered particles carry some electric charges {e i } and study scattering amplitude involving one soft photon emission with polarization and momentum (ε, k). So in this case the soft expansion of (N + 1) particle amplitude takes the following form is defined through the relation (1.2). In (1.8), the expression of tree level "soft factor" for single soft photon emission reads [5,8,[24][25][26][27][44][45][46][47][48][49][50] where the generic expression for the non-universal term N ρσ i (−p i ) contributing to the subleading soft photon theorem has been derived in [5] and its explicit form in provided in (A.3). N ρµ i (−p i ) depends on the non-minimal coupling of two finite energy fields to a soft photon field through the field strength, and the general structure of the three-point 1PI vertex involving two hard particles and a soft photon. Extending the analysis of [5] it can be argued that a complete soft factorization of order ω n for n ≥ 1 is not possible in a generic theory of QED with all possible non-minimal couplings allowed by U (1) gauge invariance, though a partial soft factorization is achievable by enforcing gauge invariance of the (N + 1)-particle amplitude [6,34]. The generalization of the tree-level soft factor S em tree for multiple soft photon emissions up to subleading order can be found in [5].
The one-loop contribution to the "soft factor" for single soft photon emission in (1.8) reads (1.11) The O(ln ω) soft factor for single soft photon emission in (1.10) has been derived in [1] as an one-loop exact result, analyzing one-loop S-matrices in the theory of minimally coupled scalar QED. There a correction to O(ln ω) soft photon factor due to gravitational interaction has also been derived. In this article we re-derive the O(ln ω) soft factor in a generic theory of quantum gravity for scattering of particles with arbitrary mass, charge and spin in presence of non-minimal coupling. This investigation will demonstrate the universal (independent of theory) nature of the O(ln ω) soft factor, while also extending the infrared divergence factorization prescription proposed in [1,51] to encompass a broad range effective field theories for charged objects. The generalization of the one-loop soft factor S em 1-loop for multiple soft photon emissions up to subleading order can be found in the section-(3.5) of [52], and the final result has been provided in (3.58).
The two-loop contribution to the "soft factor" for single soft photon emission in (1.8) reads (1.14) The O ω(ln ω) 2 soft factor for single soft photon emission in (1.12) has been derived in the section-4 of [21] as a two-loop exact result, analyzing two-loop S-matrices in the theory of minimally coupled scalar QED. From the analysis of the n-loop QED S-matrix, it is expected that the new leading non-analytic soft factor for single photon emission, as the frequency ω approaches zero, behaves like ω n−1 (ln ω) n and it relates to the tree level Nparticle amplitude. The general structure of the order ω n−1 (ln ω) n soft photon theorem is provided in [21]. The rest of the paper is organized as follows: In section-2, we establish our conventions and describe the general definition of IR-finite scattering amplitudes. We also discuss the EFT action involving massive spinning particles which transform in a generic reducible representation of the Lorentz group. In section-3, we review the covariantization prescription and define one-loop IR-finite QED S-matrices involved in the derivation of the soft photon theorem. Starting from the IR-finite S-matrices, we derive the soft photon theorem up to subleading order. In section-4, after reviewing Sen's covariantization prescription, we define the one-loop IR-finite quantum gravity S-matrices that are involved in the derivation of the soft graviton theorem. Starting from the IR-finite S-matrices, we derive the soft graviton theorem up to subleading order. At the end of both section-3 and 4, we discuss the possible generalizations of our derivations to higher orders. In section-5, we provide some open directions to explore in the future after reviewing what we have been achieved in this article.

Setup and Stretegy
Index convention: We utilize the first few Latin alphabets a, b, c, d, . . . as Lorentz indices for the tangent space, ranging from 0 to 3. The Latin alphabets starting from i, j, k, ℓ, . . . are employed as indices for identifying individual hard particles, ranging from 1 to N. The first few Greek alphabets α, β, γ, δ, . . . are used as polarization indices for spinning particles on the tangent space, while the Greek alphabets beginning with λ, µ, ν, ρ, σ, τ, . . . serve as curved space indices, ranging from 0 to 3. In section-3, where we derive the soft photon theorem solely under electromagnetic interaction in a flat background, we employ both a, b, c, d, . . . and λ, µ, ν, ρ, σ, τ, . . . as flat space Lorentz indices.
Metric and unit conventions: In our convention four dimensional Minkowski metric is η ab = diag(−1, +1, +1, +1). We work in the unit where speed of light c = 1 and Planck constant ℏ = 1 but keep the gravitational constant G explicit. We define κ ≡ √ 8πG.
Setup of scattering event: Let us consider a scattering amplitude involving N number of finite energy massive particles (hard particles) with charges, momenta, spins and polarizations {e i , p i , Σ i , ϵ i } for i = 1, 2, · · · , N and one low-energy (soft) outgoing photon/graviton with momenta and polarization k, ε, and denote this scattering amplitude by A (N +1) . In our convention we consider all the particles are incoming, so if some particles are outgoing we need to flip the sign of four momenta and electric charges for those particles. The energy of outgoing soft photon/graviton is denoted by ω so that k µ = −ωn µ where n µ being the null vector whose spatial part denotes the direction of soft photon/graviton emission. Here we are only interested to evaluate A (N +1) at one-loop order which involves Feynman diagrams involving one virtual photon/graviton running in the loop. Then we perform soft expansion (ω << |p i |) of A (N +1) to relate it with the N point amplitude which carries all the hard particles in the asymptotic state but no soft graviton, denoted by A (N ) . 4 Note that both the scattering amplitudes A (N ) and A (N +1) are distributions in momenta as A (N ) contains momentum conserving delta function δ (4) p 1 + p 2 + · · · + p N and A (N +1) contains momentum conserving delta function δ (4) p 1 + p 2 + · · · + p N + k . In four spacetime dimensions (D = 4), both scattering amplitudes exhibit infrared (IR) divergences. Therefore, our first step is to separate out the IR divergent contributions from both the scattering amplitudes in an unambiguous manner. Then we can obtain the soft factor by examining the ratio of A (N +1) and A (N ) after full/partial cancellation of the IR divergent contributions as we explain in later sections. 4 Note that the soft limit can also be defined covariantly by demanding p i .k p i .p j << 1 for all i, j = 1, · · · , N .

Feynman diagram conventions:
In all the Feynman diagrams describing scattering amplitudes, time flows from right to left and the particles involved in the scattering will always be treated ingoing. Solid lines in any diagrams corresponds to massive spinning particles and dashed lines represent photons/gravitons. If in a figure multiple Feynman diagrams appears, the counting of their numbers are always from left to right and from top to bottom. A Feynman diagram will be called an n-loop diagram only if the diagram contains n number of loops where at-least one virtual photon/graviton is propagating in each loop. The loops involving only massive virtual particles will be taken care of inside the massive EFT 1PI vertices and renormalized propagators of the massive EFT. To determine the Feynman rules for vertices involving photons/gravitons and hard particles, we follow the covariantization technique developed in the photon/graviton background in the references [2][3][4][5].
Handling IR divergences in the derivation of soft theorem: The traditional Smatrix in quantum electrodynamics and quantum gravity, in four spacetime dimensions, exhibits IR-divergence. This is due to the long-range nature of the interactions involved. Previous attempts to construct IR-finite S-matrices, beginning with the Kulish-Faddeev construction [53], demonstrated explicit cancellation of IR divergences. However, a systematic method for extracting the unambiguous IR finite part remained absent, until Grammer and Yennie provided one in [51]. A generalization of Grammer-Yennie prescription for perturbative QCD and quantum gravity can be found in [54] and [1] respectively. In the derivation of soft photon theorem, Grammer-Yennie prescription helps to factor out IR divergences from both the amplitudes A (N +1) and A (N ) in the following way (2.1) Above the exponential factor containing K em takes care of the full IR divergent contribution and the IR divergent contributions are exactly same for both the amplitudes. An explicit expression of K em is provided in (3.28). Basically Grammer-Yennie prescription provides a systematic procedure to compute IR-finite parts perturbatively for both the amplitudes. When the soft factor S em is a multiplicative function instead of differential operator, we get IR-finite . IR-finite and perform soft expansion.
In the derivation of soft graviton theorem, A (N +1) contains some extra divergent factors relative to A (N ) due to Feynman diagrams involving three graviton self-interaction vertices. An optimistic dream of the factorization IR divergence using the Grammer-Yennie decomposition proposed in [1] would be where the IR-divergent expressions of K gr and K phase are given in (4.28) and (4.29). The result mentioned above has only been verified rigorously up to one-loop order. Verifying it for all loop orders is a computationally challenging task that remains open for future investigation. Now when the soft factor S gr is a multiplicative function instead of differential operator, we get Hence to derive an unambiguous soft factor by analyzing A (N +1) IR-finite , we need to regulate the IR divergence of K phase using a cut off given by the detector resolution. This procedure can be followed to derive O(ln ω) and O ω(ln ω) 2 soft factors in (1.4) and (1.7) respectively by analyzing one and two loop IR finite amplitudes. But if we want to derive the O(ω ln ω) soft factor in the second and third lines of (1.4) we need to deal with the following additional subtleties: 1. Since the order ω ln ω soft factor in (1.4) is a differential operator we can not really commute the soft factor and infrared divergent exponential to get the second equation of (2.4). Hence to derive the O(ω ln ω) soft factor we have to start with the full divergent scattering amplitude A (N +1) instead of its IR finite part and at the end of the analysis we may be able to cancel the common IR divergent factor appears in both amplitudes in the soft theorem relation.
2. Note that the momentum conserving delta function associated with A (N +1) is On the other hand the momentum conserving delta function as- p i . Now Taylor series expansion of the first delta function around small ω produces a correction of order ω. This correction, when multiplied with the O(ln ω) soft factor, yields an additional factor of order ω ln ω at one-loop order. Therefore, this additional contribution needs to be accounted for, if it contributes something non-vanishing at this order.
In light of these additional intricacies, we have decided to postpone the derivation of the order ω ln ω soft graviton factor in (1.4) for future study and focus on deriving the order ln ω soft graviton theorem here.
EFT involving massive particles with arbitrary spin: We begin with an effective field theory (EFT) that describes the dynamics of massive spinning particles. The oneparticle irreducible (1PI) effective action for this EFT is obtained by integrating out all massive loops. The tree level amplitudes computed using this massive EFT action contain information about all the loop orders in the original un-integrated massive quantum field theory (QFT). However, if the un-integrated QFT includes massless fields, our initial approach using the 1PI effective action becomes invalid. Nevertheless, our prescription for covariantization and computation of loop amplitudes, as described below, remains valid. In such cases, the 1PI effective action should be regarded as the tree level action for the EFT. Let Φ α (x) denotes the set of all massive fields in real representation, present in the 1PI effective action 5 which transforms in a reducible representation of Lorentz group SO (1,3) in the following way, where λ ab = −λ ba is the infinitesimal Lorentz transformation parameter and Σ ab is the spin angular momentum generator of SO(1, 3) transformation in the real reducible representation. The subscript index α is used as a combined notation for denoting different fields in the theory as well as the spin/polarization indices of each of the fields. Under global U (1) EM the field Φ(x) transforms in the following way, where θ is the parameter of global U (1) EM transformation and Q is the generator of global U (1) EM transformation in the real representation of Φ(x). Usually we associate U (1) EM global charge to complex fields but since we want to covariantize the theory simultaneously in the background of gravity and gauge theory together following [5], working in terms of real field components is convenient. For example instead of a complex scalar field we work with two real scalar fields considering them in a two component vector which rotates under SO(2) and Q is the generator of SO(2) transformation. In the set of fields denoted by Φ α (x), there may be some elementary fields in the irreducible representation of Lorentz group which does not transform under global U (1) EM , for those fields the elements of the charge matrix Q will be zero. Let us start with the general form of the quadratic part of the massive particle 1PI effective action 6 S (2) = 1 2 where K(q) is the renormalized momentum space kinetic operator which satisfy the following condition: In the second lines of (2.7) and (2.8) we introduced the index free notation, which we follow through out the article. In the RHS of the above equation, + sign is for bosonic field and − sign is for fermionic field. For simplicity we work with the + sign considering Φ(x) 6 If the original theory contains some massless fields, then this action should be thought of as quadratic part of the tree level gauge fixed action. Because in presence of massless fields, the 1PI effective action of the theory may be non-local and the kinetic operator K(q) may not be polynomially expandable around q µ = 0, which is the key assumption for the validity of the covariantization prescription discussed below.
being Grassmannian even, but the final result of soft theorem computation will be same for both bosonic and fermionic fields. The Feynman propagator for the i-th particle with renormalized mass m i from (2.7) becomes where K i (q) is the kinetic term for the set of fields representing the i-th particle after proper diagonalization of the quadratic part of the action S (2) . The above equation also defines Ξ i (q) as the residue of the pole of the propagator for i-th particle. The relation between K i and Ξ i and their momentum derivatives satisfy the following relations in the index free notation, which will be useful for later computation [2] 11) 14) The Lorentz covariance of K i and Ξ i implies the following two relations where Σ i is the spin angular momentum generator for i-th component field inside Φ(x).
Taking derivatives with respect to momenta the above expressions become (2.18) Invariance of (2.7) under global U (1) EM transformation implies This also imposes constraint on the numerator of the propagator which reads The above two equations are also valid for component fields in real representation. We can take momentum derivatives on the above two relations to find useful expressions. When the i-th spinning particle is on-shell with momentum q i and polarization tensor ϵ i (q i ) it satisfies 3 Soft photon theorem at one-loop In this section we derive subleading soft photon theorem analyzing one-loop amplitudes for a quantum mechanical scattering process involving N number of massive charged particles with arbitrary spin. In [1], the order ln ω soft factor has been derived analyzing oneloop amplitude in a theory of minimally coupled scalar-QED in presence of scalar contact interaction and the soft photon factor is determined in terms of the charges and asymptotic momenta of scattered particles, and the direction cosine of soft photon emission. Here in this section we show that even for arbitrary spinning particle scattering in a generic theory of QED with non-minimal interaction, the order ln ω soft factor derived in [1] is universal (theory independent). This section should be thought of as a warm up of the next section where we are going to derive one-loop soft graviton theorem for spinning particle scattering in a generic theory of quantum gravity.

Covariantization and Feynman rules
In [5], the quadratic action S (2) in (2.7) has been covariantized simultaneously in photon and graviton background to determine 1PI vertices involving two hard spinning particles and one or two photons/gravitons up to subleading order in the expansion of the momenta of photons/gravitons. Without going into too much details, here we summarise the outcomes of the covariantization prescription in photon background and write down the Feynman rules for vertices involving one and two photons. We derive the vertices for off-shell photon with Feynman gauge fixing term, such that the Feynman propagator for virtual photon reads In position space the kinetic operator in (2.7) contains derivatives over the field Φ β (x), which have to replace by covariant derivatives under covariantization in presence of photon field A µ (x). For example in the case of one and two derivatives the covariantization rule in position space becomes Above we use the following symmetrization convention E (µ F ν) = 1 2 E µ F ν + E ν F µ for two vectors E and F . For determining minimal interaction vertices Γ (3) involving two massive spinning particles and one photon up to one derivative on photon field, and Γ (4) involving two massive spinning particles and two photons with no derivative on any of the photon fields, the information of covariantization for single and two derivatives as done above would be enough. In momentum space these covariantization rules generate the following minimally interacting actions starting from (2.7): and is the Fourier transform of gauge field defined through the following relation: At the order of one derivative on the gauge field, we can have non-minimal coupling of photon with matter fields interacting via field strength. In momentum space, the general form of non-minimal interaction takes the following form and B(q 2 ) satisfies the following relations In the last equation above + sign is for Grassmannian even field and − sign is for Grassmannian odd field. Again during the derivation we consider Φ field components being Grassmannian even, but the final result will be valid for both Grassmannian even and odd fields. We introduced index free notations for all the equations above. Starting from (3.4) and (3.6), the interaction vertex describing two ingoing spinning particles with momenta q and −(q + ℓ), polarization index α and β, and one ingoing photon with momentum ℓ becomes Expanding in small ℓ limit and using (2.8), (2.19), (3.10) the above vertex reduce to the following polarization index suppressed form Similarly starting from (3.5) the polarization index suppressed four point interaction vertex involving two incoming spinning particles with momenta q and −(q + ℓ 1 + ℓ 2 ) and two incoming photons with momenta ℓ 1 and ℓ 2 becomes We denote the scattering amplitude describing N -number of spinning hard particle scattering in massive EFT by Γ (N ) , which can be expressed as a polarization tensor contracted form in the following way where Above ϵ α i (p i ) denotes the polarization tensor for i-th spinning particle with momentum p i . Here we also should remember that Γ (N ) is a distribution as it contains a momentum conserving delta function δ (4) p 1 + p 2 + · · · + p N . We denote a part of amputated Green's function involving N number of spinning hard particles and one photon with momentum ℓ by Γ (N +1) µ (ℓ), which describes sum of the contributions of the Feynman diagrams where the photon is not connected to any external leg. Since Γ (N +1) µ (ℓ) does not contain any loop involving massless particles, we can write down the following relation between Γ (N +1) µ (ℓ) and Γ (N ) using the same covariantization prescription described above In the above expression Γ (N ) contains momentum conserving delta function δ (4) p 1 + p 2 + · · ·+p N and Γ (N +1) contains momentum conserving delta function δ (4) p 1 +p 2 +· · ·+p N +ℓ .

Grammer-Yennie decomposition and IR-finite amplitudes
In this section we discuss the Grammer-Yennie decomposition introduces in [1,51] and show how it helps to factorize the IR-divergent exponential from the IR-finite part of the amplitudes as proposed in (2.1). In Feynman gauge we decompose the photon propagator with momentum ℓ flowing from the leg i to the leg j for i ̸ = j Note that p i and p j refer to the external momenta flowing into the legs i and j, and not necessarily the momenta of the lines to which the photon propagator attaches (which may have additional contribution from external soft photon momentum or internal virtual photon momentum). For i = j we do not carry out any decomposition i.e. for the virtual photon involved in the self energy loop we do not carry out the KG-decomposition of (3.18). Also if one or both ends of the virtual photon propagator are attached to any internal massive particle propagator carrying sum of two or more external massive particle's momenta, or vertices of the massive EFT involving more than two massive particles in a loop diagram, then we do not need to perform KG-decomposition as those loop diagrams are IR-finite. In (3.18) the propagator part containing K (ij) will be denoted by K-photon propagator and the propagator part containing G (ij) will be denoted by G-photon propagator throughout this section.
Ward identities involving K-photon: Since K-photon propagator is proportional to ℓ µ ℓ ν i.e. pure gauge, we can study the Ward identity for an off-shell un-amputated three particle Green's function with one leg being the K-photon as drawn in Fig.1. The LHS of Fig.1 after contracting with ℓ µ representing K-photon can be expressed as Figure 1. This figure is a Feynman diagrammatic representation of the expression in (3.22). Solid lines represent the massive spinning particles, dashed lines represent the ingoing virtual photon with momentum ℓ and the arrow in the photon line represents that it is a K-photon (pure gauge part contracted). The solid blobs in the RHS represent a new kind of vertices and the Feynman rules for those vertices are just Q i and Q T i as written next to the vertices.
To write down the Feynman rule in the first line above for the diagram in Fig.1 i } within the square bracket above and using the property (2.13), the RHS of the above expression translates to If we un-do the small ℓ expansion in the second line of the above expression 7 and use the relation (2.20) the Ward identity turns out to be the following expression which has been diagrammatically represented in Fig.1.
Important to note that the solid blob vertices in Fig.1 carry only the information of charge of the particle with which the K-photon is interacting. The Feynman rules for the blob vertices are independent of the momenta or any other information of the theory. We also need to study the consequence of Ward identity due to insertion of a K-photon in presence of an external off-shell photon with momentum k and Lorentz index ν. The set of Feynman diagrams describing the four point un-amputated Green's function with one photon and one K-photon has been drawn in the first line of Fig.2. Using the Ward identity described in Fig.1 for first and third diagrams in the first line of Fig.2, we find the diagrams drawn after the equality in Fig.2. Now if we can show that the sum of the contribution of the three diagrams in the last line of Fig.2 vanishes, then the Ward identity of Fig.1 is also valid in presence of an external photon line. The three diagrams in the last line of Fig.2 has been drawn again in Fig.3 and the contribution becomes (3.23) Figure 2. Diagrams in the first line represents the contribution to the four point un-amputated Green's function with one photon and one K-photon. By using the Ward identity in Fig.1 we get the diagrams in second and third lines.
Substituting the expressions for the vertices from (3.12) and (3.13), the expression inside the square bracket of (3.23) turns out to be  By Taylor expanding the first line above for small ℓ and only keeping terms up to linear order in ℓ or k, we can use the identities in (2.19) and (3.8) to show that the sum of the contribution in the three lines above vanishes up to linear order in ℓ or k. This proves the diagrammatic identity in Fig.3. By utilizing the expressions of vertices from (3.12) and (3.13), which are given up to linear and zeroth order in ℓ or k respectively, it may initially appear that the validity of the result in Fig.3 is limited to linear order in ℓ or k. However, it is important to note that the results depicted in Fig.3 and Fig.1 hold true for all orders in the expansion of small ℓ and k. These results play a crucial role in establishing the gauge invariance of any amplitude involving external photons in quantum electrodynamics, as they are connected to the Ward-Takahashi identity of QED. For further details on the spinor-QED case, please look at [51] and section-(7.4) of [55]. For an unamputated Green's function with two massive spinning particles and arbitrary number of external photon legs, one insertion of K-photon in all possible way finally reduces to sum over sets of diagrams where the K-photon is connected in the end of the spinning particle legs with solid blob vertices as discussed above. This strong statement can be proved using the identity in Fig.1 and the generalized identity in Fig.4. The identity in Fig.4 is a straight forward generalization of the example discussed in Fig.3, which has been tested with the covariantized vertices up to linear order in photon momenta expansion for Γ (n) and Γ (n+1) vertices with n = 3.

IR-finite amplitudes:
As we defined earlier, A (N ) represents the all loop scattering amplitude with N number of external massive spinning particles, and A (N +1) represents the all loop scattering amplitude with N number of external massive spinning particles and one external photon. If the massive spinning particles carry definite charges {e i } then the following identity holds The K-photon insertion Ward identities of Fig.1 and Fig.4 imply the exponentiation of the one-loop K-photon contribution K em in (2.1), as proven in [51]. However, the proof of this exponentiation in [51] is valid only when the tree-level amplitude of the massive EFT Γ (N ) is independent of the momenta of scattering particles i.e. it is described by a momentaindependent contact interaction between N number of massive fields, as considered in [1]. The validity of the exponentiation of K em in a generic theory of QED, incorporating all possible interactions without considering to any kind of approximation (like assuming that virtual photon momenta are significantly smaller than external massive particle momenta), is too good to be true. 8 But in the limit of small virtual photon momenta the IR divergent piece of K em exponentiate, which is known as the leading Eikonal exponentiation. The final outcome of the Grammer-Yennie decomposition of virtual photon propagator in (3.18) is . IR-finite represent the infrared finite components of the N -particle and N -particle-1-photon amplitudes, respectively. These components are obtained by removing the exponentiated IR-divergent parts from the original divergent amplitudes defined through the relations (3.26) and (3.27) IR-finite comprise contributions from the corresponding tree-level amplitudes and loop amplitudes up to all orders in perturbation theory. However, there is a condition: if both ends of a virtual photon propagator are connected to external massive spinning particle lines (which may already contain additional real or virtual photon lines), then this photon propagator should be replaced by a G-photon propagator when we evaluate them for the IR-finite parts. Additionally the same set of diagrams need to be evaluated with Kphoton propagator as well and then have to subtract by a factor of K em times the IR finite amplitude at one less loop level. In our convention the tree level amplitudes are given by where in the subscript '0' corresponds to 0-loop i.e. tree level. At one-loop order, A IR-finite are given by: where in the subscript '1' corresponds to 1-loop. Above G,1 corresponds to the diagram in Fig.5 which we need to evaluate with G-photon propagator. non-div,1 represents the set of diagrams in Fig.11, evaluated with full photon propagator.

Derivation of soft photon theorem
The goal here will be to derive the order ω −1 and ln ω soft factors from the ratio of A IR-finite when the external photon energy is small i.e. ω << |p µ i |.

IR-finite one loop N -particle amplitude
First we want to analyze all the Feynman diagrams contributing to (3.30) and show that the result is IR-divergence free. We also provide an explicit integral expression of IR-finite 1loop amplitude for A (N ) IR-finite,1 . Using the Feynman rules derived in section-3.1, the diagram in Fig.5 with G-photon propagator contributes to the following G,1 , where the virtual photon propagator is the G-photon propagator connected between two external hard particle lines. The diagram also contributes to A (N ) K−finite,1 when evaluated with K-photon propagator and subtracted K em Γ (N ) from it's contribution.
where to get the last two lines we have used the identity in (A.9) for both i-th and j-th particles and Taylor expanded the numerator in the limit |ℓ µ | << |p µ i |, |p µ j |.
Now evaluating this expression using the identity in (A.9) and Taylor expanding Γ self,1 , where the virtual photon propagator is the full photon propagator. In the last diagram the cross corresponds to counter term, which cancels the UV divergences. non-div,1 , where the virtual photon propagator is the full photon propagator whose one or both end connected to Γ (N ) . In the last diagram the cross corresponds to counter term, which cancels the UV divergences. Fig.6 renormalizes the massive spinning particle propagators in presence of electromagnetic interaction, and all the loops are IR-finite. Say the three diagrams in Fig.6 contributes to

Diagrams in
respectively. Where F 1 , F 2 , C are unknown constant matrices, which are related using onshell renormalization condition Hence the on-shell renormalization condition implies In Fig.7 we draw the sets of diagrams where one or both ends of the photon loop are attached to some internal massive virtual line or massive EFT vertex inside Γ (N ) . These diagrams are also IR-finite. For example the first diagram in Fig.7 with full photon propagator becomes Now using the identity in (A.9) for i-th particle and the Feynman rule of (3.17) it is evident that in the limit when the loop momentum ℓ µ → 0 the integration behaves as d 4 ℓ |ℓ| 3 , hence is IR-finite. This is the reason we call the set of diagrams in Fig.7 as A (N ) non-div,1 as those are IR-divergence free. We do not need to evaluate the contribution A (N ) non-div,1 explicitly for deriving soft photon theorem. Now summing over the contributions of (3.32), (3.34), (3.37) and A (N ) From the above expression, it is clear that the loop integrals become infrared finite in the limit as ℓ approaches zero. Therefore, the Grammer-Yennie prescription offers a clear definition of an IR-finite S-matrix. An explicit expression of IR-finite S-matrix at one-loop order is provided in the above expression in a general theory of quantum electrodynamics (QED).

IR-finite one loop (N + 1)-particle amplitude in the soft limit
Here we analyze all the Feynman diagrams contributing to (3.31) in the soft limit i.e. ω → 0. Let us start analyzing the first diagram in Fig.8 with G-photon propagator, which has the following expression after using Feynman rules Substituting the results from (A.4) and (A.9) in the above expression and after some manipulation we get , where the virtual photon propagator is the G-photon propagator connected between two external particle lines. Here we omitted the diagrams involving counter terms to remove UV divergences. We need to sum over all possible external legs while evaluating the contributions from these diagrams. When we evaluate these diagrams with K-photon propagator and subtracted K em Γ (N +1) from it's contribution, it also contributes to A (N +1) Inside the square bracket of the numerator in the above expression, we only keep the terms up to linear order in ℓ or k, as our vertices are derived only up to that order. Note that the above expression is IR-finite in the limit ℓ µ → 0 and in the region of integration |ℓ µ | << ω << |p µ i |, |p µ j | it contributes at order O(ω 0 ). Now to extract ln ω contribution we approximate the integrand in the integration range ω << |ℓ µ | << |p µ i |, |p µ j |, by doing so we can approximate With this approximation the order ω −1 and ln ω contribution turns out to be, Above the subscript "reg" in the loop integral corresponds to the restricted loop-momentum range ω << |ℓ µ | << |p µ i |, |p µ j |. Using Feynman rules the second diagram in Fig.8 with G-photon propagator becomes After substituting the results from (A.6) and (A.9) and keeping terms which can contribute up to order ln ω in the integration range reg ≡ (ω << |ℓ µ | << |p µ i |, |p µ j |) we get, self,1 , where the virtual photon propagator is the full photon propagator connecting two different points on the same massive spinning particle leg. The cross appears in some diagrams above corresponds to counter term, which cancels the UV divergences in the renormalization prescription.
Using Feynman rules the third diagram in Fig.8 with G-photon propagator becomes After substituting the results from (A.7) and (A.9) and keeping terms which can contribute up to order ln ω in the integration range ω << |ℓ µ | << |p µ i |, |p µ j | we get, The fourth diagram in Fig.8 with G-graviton contributes at order ω −1 but won't contribute at order ln ω. The order ω −1 contribution turns out to be, The fifth diagram in Fig.8 start contributing from order ω 0 in the soft expansion when we evaluate it with G-photon propagator i.e. A V = 0 + O(ω 0 ). Here we are not writing down the non-vanishing contribution of A V at order ω 0 explicitly, as it is not essential for deriving the order ω −1 and ln ω soft factors. Now summing over the external particle legs, total contribution of A (N +1) G,1 at orders ω −1 and ln ω turns out to be Interestingly in the above expression all the theory dependent pieces involving K i , Ξ i , B i at order ln ω cancels out at the integrand level, when we some over the contributions (3.43),(3.45), (3.47). This confirms the fact that the subleading soft photon theorem at order ln ω is universal. The above result can also be rewritten in the following compact form 50) Figure 10. The first diagram above corresponds to the sum over the contributions of the first four diagrams in Fig.9 when we replace the photon polarization ε µ by k µ . The second diagram above corresponds to the sum over the contributions of the fifth and sixth diagrams in Fig.9 when we replace the photon polarization ε µ by k µ .
where the expression of A (N ) G,1 is given in (3.32), and K reg em is the approximated form of the integral K em in (3.28) in the integration range ω << |ℓ µ | << |p µ i |, |p µ j |. The integration has been explicitly evaluated in [1] and the result reads (3.51) Above η j = +1 if j-th particle is ingoing and η j = −1 if j-th particle is outgoing and under ≃ sign we only have written the order ln ω contribution while evaluating the integral.
To evaluate the contribution of A where the expression of A (N ) K−finite,1 is given in (3.34). Note that the above result contributes at order ω −1 and does not contribute to order ln ω soft theorem.
The self-energy kind of diagrams contributing to A (N +1) self,1 in Fig.9 are not necessary to compute explicitly as these diagrams sum up to zero using on-shell renormalization condition (3.36) as we are going to discuss below. The contribution from the sum of the first four diagrams in Fig.9 can be described by the following general structure ϵ T i ε.
. Similarly the contribution from the sum of fifth and sixth diagrams in Fig.9 can be described by the following structure ϵ T i ε.p i f 2 (p i .k)Γ (N ) (i) (p i + k). Where f 1 (p i .k) and f 2 (p i .k) are two unknown functions with specific polarization/spin indices, which we determine below by replacing ε µ → k µ and using Ward identity. Using the diagrammatic identities of Fig.1 and Fig.3, the first four diagrams in Fig.9 after replacing ε µ → k µ reduces to the first diagram in Fig.10. Similarly using the diagrammatic identities of Fig.1 and Fig.4, the sum of fifth and sixth diagrams in Fig.9 after replacing ε µ → k µ reduces to the second diagram in Fig.10. Now using the constant matrices F 1 , F 2 introduced in (3.35) for the diagrams in Fig.6 and comparing the general structure above we get Now we substituting the above result in the general structures for the sum of diagrams mentioned above. Finally summing over the contributions of the first seven diagrams in Fig.9 we get Note that F 1 , F 2 satisfy the same property under the operation of charge matrix Q as Ξ i satisfies in (2.20). Hence using this relation the terms inside the square bracket in the above expression reduces to Q T i (F 1 + F 2 + C) which vanishes using the on-shell renormalization condition (3.36). Note that in the above expression we neglected the possible order ω 0 contribution in the soft expansion of counter term diagram. Using the same on-shell renormalization condition the sum of the rest of the diagrams in Fig.9 also vanishes up to possible order ω 0 contribution . Hence in the soft limit, the sum over all the Feynman diagrams in Fig.9 contribute as (3.55) The diagrams in Fig.11 contributing to A (N +1) non-div,1 are IR-divergence free, as in ℓ µ → 0 limit and finite k µ the third, sixth and seventh diagrams behaves like d 4 ℓ |ℓ| 3 , and first, second, fourth and fifth diagrams behaves like d 4 ℓ |ℓ| 2 . On the other hand the sum of first, fourth and sixth diagrams in Fig.11 contribute to leading soft factor at order ω −1 . Now in the integration region ω << |ℓ µ | << |p µ i |, individually the first and third diagrams in Fig.11 behave like reg d 4 ℓ ℓ 2 −iϵ 1 (p i .ℓ+iϵ) 2 after the expansion of the propagators, hence those have the possibility of contributing at order ln ω. But when we sum over the contributions of the first, second and third diagrams, the order ln ω contributing coefficient of reg cancels each other, and left out part starts contributing from order ω 0 . Hence leaving first, fourth and sixth diagrams, all the other diagrams start contributing at order ω 0 in the soft expansion. After summing over all the contributions, we get non-div,1 + O(ω 0 ) . (3.56) Figure 11.
Set of 1-loop diagrams contributing to A (N +1) non-div,1 consists of diagrams where the virtual photon propagator is the full photon propagator, with at least one leg connected to an internal massive particle propagator or massive EFT vertices. Diagrams involving counter terms to remove UV divergences have been omitted.
Soft photon theorem result at one-loop: Summing over the contributions of (3.50),(3.52),(3.55),(3.56) in the soft limit, we get the following soft theorem expression This result agrees with the loop corrected subleading soft photon theorem, originally derived in [1] for minimally coupled scalar QED. The derivation of this result further confirms the universality of the ln ω soft factor in scattering events involving particles with arbitrary spins in a generic theory of quantum electrodynamics, which allows for arbitrary non-minimal couplings. Moreover, the obtained result also verifies the well-known fact that Weinberg's leading soft photon theorem remains unaltered by loop corrections. Therefore, even from the analysis presented above, we observe that Weinberg's soft theorem still holds, relating two one-loop IR-finite amplitudes.

Discussion on generalization
In order to obtain the order ω ln ω soft factor from the aforementioned one-loop amplitude, we require the vertices: Γ (3) as given in equation , new sets of non-minimal couplings will contribute to Γ (3) . Due to these reasons, we are unable to derive the order ω ln ω soft photon theorem in this article, and it is not clear whether such a soft factorization at order ω ln ω is possible or not.
In the theory of scalar-QED minimally coupled to gravity the order ln ω correction to (3.57) due to gravitational interaction has also been derived in [1] by analyzing one-loop amplitudes. In the theory of scalar-QED soft photon theorem at order ω(ln ω) 2 has also been derived in [21] analyzing two-loop amplitudes and the soft factor is provided in (1.12). Multiple soft photon theorem up to subleading order in soft expansion has also been derived in section-(3.5) of [52], by analyzing one-loop amplitudes and the result reads:

Soft graviton theorem at one-loop
In this section we derive subleading soft graviton theorem analyzing one-loop amplitudes for a quantum mechanical scattering process involving N number of massive particles with arbitrary spin and one graviton, extending the analysis of [1] for a generic theory of quantum gravity. This derivation will establish the universal (theory independent) nature of ln ω soft graviton factor.

Sen's Covariantization prescription and Feynman rules
The covariantization of the quadratic part of the massive EFT action (2.7) in the soft gravitational background has been carried out in [2][3][4]. This development has been utilized to derive vertices involving two massive spinning particles and one or two on-shell soft gravitons. However, when performing loop computations, we require the same vertices involving off-shell gravitons. Due to the off-shell nature of the gravitons, it is not possible to independently impose the traceless and transverse conditions on the gravitational fluctuation consistently with diffeomorphism. To address this issue, we make slight modifications to the covariantization prescription proposed in [2][3][4], as described below. In this work, we use a different parametrization of gravitational fluctuation, which is closely related to the one presented in [56]. This alternative parametrization allows us to derive the vertices within the covariantization prescription while ensuring compatibility with the de Donder gauge choice. Let us define the deviation of background metric from flat Minkowski metric as where κ = √ 8πG with G being the four dimensional Newton's constant. In de Donder gauge ∂ µ h µν = 1 2 ∂ ν h ρ ρ the graviton propagator becomes where (µν) and (ρσ) are Lorentz indices of the two ends of the graviton propagator. The polarization tensor for on-shell graviton with momentum k will be denoted by ε µν (k) which satisfies the traceless and transverse conditions With the definition of gravitational fluctuation (4.1), we express the following quantities as a power series expansion of κ: In the RHS of all the above expressions both curved space indices (µ, ν, ρ, σ, . . .) and tangent space indices (a, b, c, . . .) are raised or lowered by using the Minkowski metric η. The trace of the metric fluctuation is defined as h ≡ h µν η µν . Additionally, the symbols e a µ represent the vierbein, E µ a represents the inverse vierbein, ω ab µ represents the spin-connection, Γ λ µν represents the Christoffel connection, and R µνρσ represents the Riemann tensor.
Covariantization: In the covariantization prescription, we derive the action that describes the interaction between two spinning particles and one graviton up to second derivative order on gravitational fluctuations. Additionally, we obtain the interaction between two spinning particles and two gravitons up to first derivative order on gravitational fluctuations in the derivative expansion. At the second derivative order, the interacting action of two spinning particles and one graviton also includes a generic non-minimal coupling through the curvature tensor, which is inherently general coordinate invariant. Under covariantization prescription tangent space derivatives on Φ β should be replaced by covariant derivative in the curved space after multiplication of inverse vierbeins in the following way where the expression of one, two and three covariant derivatives on Φ β are given by The last equation above needs to symmetrize in µ, ν, ρ indices before substituting in (4.5).
In the above expressions the terms within square brackets are new at each derivative order in the covariantization prescription, which are important in deriving the interaction between two spinning particles and one graviton up to two derivative order on gravitational fluctuation, and the interaction between two spinning particles and two gravitons up to one derivative order on gravitational fluctuation. To the derivative order we are interested to find the interacting parts of the action we do not need to know the new terms coming from covariantization of more than three derivatives on Φ. Also while covariantizing the action (2.7), we need to include √ −det g as a covariant measure of volume. The Fourier transform of gravitational fluctuation h µν (x) will be denoted by h µν (ℓ) and the relation between them is given by (4.9) Interaction part of the action: Under the above prescribed covariantization procedure the interactive part of the action describing interaction between two massive spinning particles and one graviton up to quadratic order in graviton momentum becomes By covariantization procedure we only get the interacting action representing minimal coupling of graviton with matter field. On top of it at the quadratic order in graviton momentum we also need to add generic non-minimal interaction term which describes the interaction between two spinning particles and one graviton through lineariszed Riemann tensor. A generic form of the action describing such kind of non-minimal interaction is given by where R µνρσ (ℓ) is the Fourier transform of the linearized Riemann tensor in (4.4) which reads In the non-minimal action (4.11), G satisfies the following property G αβ,µνρσ (q 2 ) = ± G βα,µνρσ (−q 1 − ℓ) . (4.13) In the above equation + sign is for Grassmannian even field and − sign is for Grassmannian odd field. Again during the derivation we consider Φ field components being Grassmannian even, but the final result will be valid for both Grassmannian even and odd fields. Following the covariantization procedure, the part of the action describing interaction between two massive spinning particles and two gravitons up to linear order in graviton momenta becomes (4.14) We also need to provide a purely gravitational effective action, constructed out off curvature tensors in derivative expansion. This action describes the self-interaction of the graviton field and provide dynamics to graviton. For the analysis we are conducting here, it suffices to consider the leading term of the EFT action, which corresponds to the Einstein-Hilbert action and is expressed as follows: Feynman rules for the vertices: Starting from the interacting parts of the action given in (4.10) and (4.11), Feynman rule for the vertex describing interaction between two spinning particles with momenta q and −(q + ℓ), and one graviton with Lorentz indices µν and momentum ℓ turns out to be where we suppressed the massive particle spin/polarization indices. The above expression is symmetrized under µ ↔ ν exchange and in our convention momenta of the particles are always flowing towards the interaction vertex i.e. ingoing. Analogously starting from the interacting part of the action given in (4.14), Feynman rule for the vertex describing interaction between two spinning particles with momenta q and −(q + ℓ 1 + ℓ 2 ), and two gravitons with Lorentz indices (µν), (ρσ) and momenta ℓ 1 , ℓ 2 respectively, turns out to be where we suppressed the massive particle spin/polarization indices. We could have symmetrize the above expression under µ ↔ ν and ρ ↔ σ exchanges. However, it is unnecessary because any contraction involved with this vertex in the computation of loop diagrams will already exhibit symmetry under these exchanges.

KG-decomposition and IR-finite amplitudes
In Grammer-Yennie prescription we decompose the internal graviton propagator (4.2) with momentum ℓ flowing from massive spinning particle i to j(̸ = i) in the following way [1] (see also [57]) Note that p i and p j above refer to the external momenta flowing into the legs i and j, and not necessarily the momenta of the lines to which the graviton propagator attaches (which may have additional contribution from external soft graviton momentum or internal virtual graviton momentum). For virtual gravitons whose one or both ends are attached to a 3-graviton vertex instead of a massive particle, or to some internal massive particle line or vertex inside Γ (N +1) , we do not carry out any Grammer-Yennie decomposition as those won't contribute to IR divergences. In (4.20) the propagator part containing K (ij) will be denoted by K-graviton propagator and the propagator part containing G (ij) will be denoted by G-graviton propagator throughout this section.
Ward identities involving K-graviton: From the definition of K-graviton propagator in (4.21), it is clear that K-graviton is proportional to a pure gauge of structure ζ µ ℓ ν + ζ ν ℓ µ , with ζ = p i −ℓ when it flows from i-th leg. Let us study the Ward identity for an off-shell unamputated three particle Green's function involving two massive fields and one K-graviton. It is diagrammatically represented in Fig.12. Figure 12. This figure is a Feynman diagrammatic representation of the expression in (4.25). Solid lines represent the massive spinning particles, dashed lines represent the ingoing virtual graviton with momentum ℓ and the arrow in the graviton line represents that it is a K-graviton proportional to ζ µ ℓ ν + ζ ν ℓ µ . The solid blobs in the RHS represent a new kind of vertices and the Feynman rules for the left blob vertex is −κ{2ζ.(q +ℓ)+ζ ν ℓ b Σ νb } and for the right blob vertex is κ{2ζ.q −ζ ν ℓ b Σ T νb }.
The LHS of the Fig.12 takes the following form: To compute the above expression we used the result of (C.4) and then simplified using the identities in (2.10),(2.17) and derivatives of (2.17). After all the simplification we get Now undoing the ℓ expansion for Ξ(−q − ℓ) the above expression can be re-written as (4.25) The above identity is diagrammatically represented in Fig.12. It is worth noting that the blob vertices in this representation depend on the momenta of the massive particle and the attached K-graviton, as well as on the spin angular momenta of the massive particle, as indicated in Fig.12.
In comparison to the QED Ward identity in (3.22), the momenta dependence of the blob vertices weakens the power of KG-decomposition. For instance, if we wish to study the Ward identity for the four-point un-amputated Green's function associated with two massive particles, one graviton, and one K-graviton, analogous to Fig.2, it will result in new Feynman rules for the right blob vertex drawn in the second diagram after the equality of Fig.2. Specifically, the right blob vertex Feynman rule reads κ{2ζ.(q + k) − ζ ν ℓ b Σ T νb }, which contains an additional term 2κζ · k relative to the three-point un-amputated Green's function right blob Feynman rule κ{2ζ.q − ζ ν ℓ b Σ T νb } in Fig.12. In turn this implies that A (N +1) contains an extra exponentiation term compared to the IR-exponentiation factor of A (N ) . However this extra contribution arising from the momenta dependent right blob vertex rule contributes to IR-finite part in the loop integral, allowing us to follow the same strategy used in deriving the soft photon theorem. In [1], it is also observed that the righthand side of Fig.3 in presence of one external graviton, due to one K-graviton insertion, does not vanish, but instead leaves out some finite residual contribution. Fortunately, this residual terms also contributes to the IR-finite part in the loop integral when we evaluate A (N +1) , which we can systematically account for it as well. However, as mentioned in section-2, we still need to regulate some additional IR divergence in A (N +1) that arises from the Feynman diagram involving three graviton interaction vertices of which one graviton being the external graviton at the one-loop level.

IR-finite amplitudes:
The IR-finite amplitudes associated with the scattering of N number of spinning massive particles, and associated with N number of spinning massive particles plus one outgoing graviton are defined by is the Eikonal IR-divergent exponentiated factor. In equations (4.26) and ( IR-finite represent the infrared finite components of the N -particle and N -particle-1-graviton amplitudes, respectively. 9 These components are obtained by removing the exponentiated IR-divergent parts from the original divergent amplitudes defined through the relations (4.26) and (4.27). Both A IR-finite , we need to use an explicit IR cut-off for the diagrams involving graviton self interaction vertices. The "reg" over the ≡ sign in (4.27) corresponds to this particular IR regularization scheme. Effectively this IR-regularization scheme removes a factor of exp{K phase } from A (N +1) with an explicit expression for K phase being corresponding tree-level amplitudes and loop amplitudes up to all orders in perturbation theory. However, there is a condition: if both ends of a virtual graviton propagator are attached to different external massive spinning particle lines (which may already contain additional real or virtual graviton lines), then this graviton propagator should be replaced by a G-graviton propagator when we evaluate them for the IR-finite parts. Additionally the same set of diagrams need to be evaluated with K-graviton propagator as well and then have to subtract by a factor of K gr times the IR finite amplitude at one less loop level. On top of it to evaluate the full IR-finite part of A (N +1) , we need to IR regulate diagrams containing at least one graviton self interacting vertex involving the external graviton. Rigorous definitions of one-loop IR-finite amplitudes are provided below. Now analogous to (3.29), (3.30) and (3.31), here we also split the IR-finite parts of the tree and one-loop amplitudes in the following way non-div,1 .

A (N )
non-div,1 represents the set of diagrams in Fig.7, evaluated with full graviton propagator representing the dashed lines.

A
(N +1) G,1 corresponds to the set of diagrams in Fig.8 where we need to evaluate the diagrams with G-graviton propagator representing the dashed virtual lines. The dashed external lines represent the on-shell graviton with momentum k.

A (N +1)
K−finite,1 corresponds to the contribution from the diagrams in Fig.8, evaluated with K-graviton propagator representing the dashed line and then subtracted the contribution K gr × Γ (N +1) at the integrand level. The dashed external lines represent the on-shell graviton with momentum k.

A (N +1)
3-graviton-reg,1 represents the sum of the contributions of Feynman diagrams in Fig.13, when we evaluate them using full graviton propagator and regulate the IR divergence considering detector resolution as the explicit IR cut-off.

A
(N +1) self,1 represents the set of Feynman diagrams in Fig.14, evaluated with full graviton propagator representing the dashed lines. The dashed external lines represent the onshell graviton with momentum k.

A
(N +1) non-div,1 represents the set of diagrams in Fig.11, evaluated with full graviton propagator representing the dashed lines. The dashed external lines represent the on-shell graviton with momentum k.

Derivation of soft graviton theorem
The goal here will be to derive the order ω −1 and ln ω soft factors from the ratio of A IR-finite when the external graviton energy is small i.e. ω << |p µ i |.

IR-finite one-loop N-particle amplitude
Let us evaluate the sum of contributions from Fig.5 with insertion of G-graviton propagator and the finite part with the insertion of K-graviton propagator together. We can do this by examining the diagram in Fig.5, where the dashed line represents a full graviton propagator. Then we need to subtract K gr × Γ (N ) from the evaluated result. The mathematical expression reads (4.34) We evaluate the above expression by using the identity derived in (C.4). The result in small ℓ expansion in the integrand turns out to be Note that in the limit ℓ µ → 0 the integrand of the above expression behaves like d 4 ℓ |ℓ| 3 at leading order, hence the contribution is IR-finite as promised. Fig.6 renormalizes the massive spinning particle propagators in presence of gravitational interaction, when the dashed lines represent graviton propagators. The contribution from the loop diagrams there are IR-finite. Hence following the analogous wave function renormalization condition (3.36), the sum of the contribution vanishes i.e.

Diagrams in
Now let us analyze the diagrams in Fig.7 with dashed lines being full graviton propagators are connected to some internal massive virtual lines or massive EFT vertices inside Γ (N ) . In the limit ℓ µ → 0 the integrand of the first diagram behaves like d 4 ℓ |ℓ| 3 at leading order, and the second diagram behaves like d 4 ℓ |ℓ| 2 at leading order, hence IR-finite. We do not need to evaluate them explicitly. Let the total contribution after removing the UV divergences by adding counter term diagrams reads non-div,1 . (4.37) Hence the total IR-finite contribution to N -particle amplitude follows from the definition (4.32) becomes non-div,1 . (4.38)

IR-finite one-loop (N+1)-particle amplitude in the soft limit
Here we analyze all the Feynman diagrams contributing to (3.31) in the soft limit i.e. ω → 0. We start by analyzing the diagrams in Fig.8 with the dashed lines being full graviton propagators, which evaluates the contribution of A K−finite,1 +K gr Γ (N +1) . By evaluating this sum with full graviton propagator, we avoid all the computational complicacies in the KG-decomposition in presence of external graviton as pointed out in the paragraph below (4.25). Then finally we subtract the contribution of K gr Γ (N +1) from the sum to extract the ω −1 and ln ω soft factors from A K−finite,1 in the limit ω << |p µ i |. The first diagram in Fig.8 with full graviton propagator representing the internal dashed line takes the following form After using the identities in (C.4) and (C.6) and simplifying the above expression reduces to The second diagram in Fig.8 with full graviton propagator representing the dashed internal line takes the following form After using the identities in (C.4) and (C.8) and simplifying the above expression reduces to The third diagram in Fig.8 with full graviton propagator representing the dashed internal line takes the following form After using the identities in (C.4) and (C.6) and simplifying the above expression reduces to The fourth diagram in Fig.8 with full graviton propagator representing the dashed internal line takes the following form 45) The fifth diagram in Fig.8 with full graviton propagator representing the dashed internal line takes the following form In both the expressions of B IV and B V the integrands can be simplified using (C.4) and the result of the following common loop integral reads Hence the total IR-finite contribution from the sets of diagrams in Fig.8 reads where the (N + 1) particle tree level amplitude is given by (4.49) In the soft limit the above tree level (N + 1) particle amplitude provides the tree level soft graviton theorem with soft factor given in (1.3). Substituting the results of (4.40), (4.42), (4.44), (4.45), (4.46) and (4.28) in (4.48), it is easy to see that the final expression is IR-finite in the limit ℓ µ → 0. In explicit computation, the IR divergent contribution in the sum of (B I + B III ) + B IV + B V cancels with the IR divergent contribution of K gr × Γ (N +1) , and B II is IR finite at finite value of k. Now let us analyze the expression (4.48) in the soft limit i.e. ω << |p µ i |. In this limit a part of the final expression (4.48) contributes to ω −1 for the full integration range of virtual graviton momenta, which reproduce the Weinberg's soft graviton factor relating two IR-finite one-loop amplitudes. For evaluating the rest of the IR-finite part of (4.48) we divide the virtual momenta integration range into three regions: |ℓ µ | ∈ [0, ω], [ω, |p µ i |] and [|p µ i |, ∞). It turns out that the integrand starts contributing at order ω 0 in the region of integration |ℓ µ | ∈ [0, ω] and [|p µ i |, ∞) 10 in the soft limit. Finally in the region of integration |ℓ µ | ∈ [ω, |p µ i |] the integrand contributes at order ln ω which is dominant compare to order 10 In the integration region |ℓ µ | ∈ [|p µ i |, ∞) the integrand is UV divergent, so needs to use UV regulator and add appropriate counter terms to extract finite result. ω 0 in the soft expansion. Hence the order ω −1 and ln ω contribution from (4.48) turns out to be Above the expression of A (N ) K−finite,1 is given in (4.35). The "reg" in the subscript of loop momentum integration above refers to the integration range |ℓ µ | ∈ [ω, |p µ i |], |p µ j |]. It is intriguing to observe that the components reliant on theory, such as K i , Ξ i , and ∂Γ (N ) ∂p i , along with the dependence on spin angular momenta of the massive particles, vanish when the individual diagram contributions are summed in the ln ω order soft factor contributing integrand. This theory and spin independence feature is not true for the order ω 0 contribution from (4.48), which we ignored here.
Let us proceed to calculate the contribution of A (N +1) 3-graviton-reg,1 . This term represents the sum of contributions from the Feynman diagrams shown in Fig.13. We evaluate them using full graviton propagator and regulate the IR divergence of the virtual loop momentum integration by introducing an explicit IR cut-off R −1 , which serves as the resolution of the detector. The expression of the first diagram in Fig.13 reads Now using the identity in (C.4) and substituting the 3-graviton vertex from (4.19) with on-shell and transverse-traceless condition for the external graviton the above expression 3-graviton-reg,1 after regulating the IR-divergence considering detector resolution as the explicit IR cut-off. The solid lines represent massive spinning particles and the dashed lines represent gravitons.
In the region |ℓ µ | ∈ [R −1 , ω] we approximate the propagator denominators at leading order as Then only when we choose the order O(kk) terms from the V (3) vertex, after the loop momentum integration we can have O(ln ω) contribution. On the other hand in the region "reg"≡ |ℓ µ | ∈ [ω, |p µ i |] we approximate the propagator denominators at leading order as Then only when we choose the order O(ℓℓ) terms from the V (3) vertex, after the loop momentum integration we can have O(ln ω) contribution. Hence the full contribution at order ln ω from all the three regions of integration becomes Note that the integrands above are independent of the spin angular momenta of external massive particles as well as do not depend on the theory dependent terms such as K i , Ξ i or non-minimal couplings. The first integrand above can be evaluated using the result of the integral (derived in [1]) Above η i convention is the same as described below (1.6). We won't evaluate the second integrand in (4.55) explicitly at this moment, but will simplify some of the terms which contains (p i .ℓ) or (p j .ℓ) in the numerator using the momentum conservation relation in (4.57). Also the term containing (ℓ.ε.ℓ) can be simplified by the following using integration by parts relation Following all the steps outlined above and using exchange symmetry (p i , ℓ) ↔ (p j , −ℓ) in the second integrand of (4.55), the simplified expression of C I becomes Note that in the first line above j = i sum is included while it was not present in the expression (4.55). The inclusion of this term originates from the second term in the numerator of the first integral in (4.55) after using the momentum conservation relation (4.57).
For completeness here we also briefly analyze the second diagram in Fig.13 which reads (4.60) In the integration region |ℓ µ | ∈ [R −1 , ω] or |(ℓ + k) µ | ∈ [R −1 , ω] the numerator of the potentially contributing ln ω terms vanishes in the integrand. On the other hand, in the integration region "reg" using a set of integration by parts to cancel (p i .k) −1 factor, the terms potentially contributing at order ln ω becomes The above expression of C II can be simplified using the following identity (4.62) and the simplified expression reads The third, fourth and fifth diagrams in Fig.13 do not contribute to order ln ω from any integration region. Hence the total order ln ω contribution to A (N +1) 3-graviton-reg,1 after summing over (4.59) and (4.63) turns out to be Let us now analyze the set of Feynman diagrams in Fig.14 which contributes to A (N +1) self,1 . All the diagrams are IR-finite for finite k and the UV divergences in the sum of contributions cancel by using on-shell renormalization condition with proper choice of counter terms. After renormalization a finite contribution remains, and it contributes at order ln ω in the region "reg". In specific, the sum of the diagrams in Set 1, Set 2, Set 4 and Set 5 individually vanishes using the gravitational analogue of the wave functional renormalization condition (3.36). On the other hand the counter term choice of the last diagram in Set 3 cancels the sum of UV divergences appearing in the the first four diagrams in Set 3 and the second and third diagrams of Fig.13. After the cancellation of these UV divergences, only in the integration region "reg" we get the following expression, which can potentially contribute at order ln ω in the soft limit Before even evaluating A (N +1) non-div,1 from the set of diagrams in Fig.11, it becomes evident that these diagrams are IR-finite for finite k when ℓ µ → 0. The sum of first, fourth and self,1 , where the virtual dashed lines represent the full graviton propagator connecting two different points on the same massive spinning particle leg. The cross appears in some diagrams above corresponds to counter term, which cancels the UV divergences in the renormalization prescription.
sixth diagrams in Fig.11 contribute to leading soft graviton factor at order ω −1 multiplying A (N ) non-div,1 . On the other hand in the integration region ω << |ℓ µ | << |p µ i |, though individually the first and third diagrams in Fig.11 have the potential to contribute at order ln ω as those behave like ω d 4 ℓ |ℓ| 4 , when we sum them up, in the numerator we get an extra factors of ℓ and/or k. Hence they can only contribute from order ω 0 or ω ln ω. Hence, summing over all the contribution we get Summing over the contribution of (4.50), (4.64), (4.65) and (4.66) we get where the expression of A (N ) IR-finite,1 is given in (4.38).
Soft graviton theorem result at one-loop: The expression for one-loop amplitude in the soft expansion derived in (4.67) can be re-written in the following compact way where the expressions of K reg phase and K reg gr are given in (1.6) and (1.5) respectively. This result agrees with the loop corrected subleading soft graviton theorem, originally derived in [1] for minimally coupled scalar-gravity with scalar contact interaction. Here the rederivation of this result from a scattering amplitude involving particles with arbitrary spins in a generic theory of quantum gravity confirms the universality of the ln ω soft factor. Moreover, the obtained result also verifies the well-known fact that Weinberg's leading soft graviton theorem remains unaltered by loop corrections. Therefore, even from the analysis presented above, we observe that Weinberg's soft theorem still holds, relating two one-loop IR-finite amplitudes.

Discussion on generalization
When some of the massive spinning particles carry electric charge, the order ln ω soft graviton factor in (4.68) undergoes correction due to electromagnetic interactions. The correction term has been derived in [1] in a minimally coupled charged scalar theory. A straightforward generalization of our derivation, combined with the covariantization prescription outlined in [5], will be useful for establishing the universal nature of this correction. At one-loop order, it is also possible to derive the spin-dependent order ω ln ω soft graviton factor as provided in (1.4), once we understand how to account for the additional subtleties highlighted in the discussion below (2.4). Furthermore, through an analysis of two-loop amplitudes in the soft limit, it is possible to derive the order ω(ln ω) 2 soft graviton theorem as conjectured in (1.7), which is also expected to be universal. We intend to pursue these analyses in the future.
Using the Feynman rules derived in section-4.1, in combination with the identities in appendix-C, it is possible to derive the tree-level simultaneous and consecutive double softgraviton theorems up to sub-subleading orders in a generic theory of quantum gravity. This result does not exist in literature and will be interesting to explore in future.

Summary and outlook
In this article, we have proved that the order ln ω soft photon and graviton theorems are universal (theory independent) by working within a generic setup. Along the way, we have also provided the definitions of IR-finite amplitudes in the generic theory of QED and quantum gravity for scattering involving spinning particles. We used the Ward Identities along with the Grammer-Yennie decomposition to perform Eikonal exponentialization and extract the IR-finite amplitudes. At one-loop order, we provide integral expressions for the IR-finite amplitudes, which have been used to derive soft factors. While it would have been desirable to explicitly evaluate these one-loop IR-finite amplitudes and explore their crossing and unitarity properties in a specific theory, this remains an avenue for future investigation. 11 We have also discussed that while the KG-decomposition is powerful for extracting IR-finite amplitudes in a generic theory of QED, it loses efficacy in the generic theory of quantum gravity. Furthermore, we have provided a set of soft photon and graviton theorems up to two-loop orders in section-1, with some of these theorems conjectured based on classical analysis.
In [13,14], a relation between low-frequency electromagnetic/gravitational waveforms and the classical limit of soft photon/graviton factors has been derived. In frequency space, the low-frequency gravitational waveform in D = 4 is given by 12 where S gr classical represents the classical limit of quantum soft factor for single soft graviton emission with momentum k µ = −ωn µ . In the above expression e µν is defined by Note that to derive the classical limit of the quantum "soft factor" from (1.1), first we need to incorporate the sum over particles indexed by i. Subsequently, in their respective expressions, when the orbital momentum operator operates on A (N ) , it must be substituted with the classical orbital momenta of the scattered objects. Furthermore, based on explicit 11 BS acknowledges Waël Aoun for pursuing his master thesis project on this topic using the IR-finite S-matrix formalism proposed in [58] and making substantial progress. 12 The explicit expression is given in equation (2.5) of [23], where we have adopted the convention that the energy of the outgoing particle is positive, which is contrary to the convention used in this paper. After Fourier transformation in frequency variable, it provides gravitational memory along with multiple tail contributions at late and early time for generic gravitation scattering event.
classical computations of gravitational waveforms presented in [16,21,23], it has been conjectured that the validity of the relation (5.1) holds true solely when, in the classical limit, the soft graviton factor is calculated using the retarded propagator for the graviton field instead of Feynman propagator. In practice this prescription suggests that in the integral representations of (1.5) and (1.6), the term (ℓ 2 − iϵ) −1 should be substituted with − (ℓ 0 + iϵ) 2 − ⃗ ℓ 2 −1 in order to extract S gr classical . With this substitution, S gr classical only receives contributions from the part of the loop integrals in which the virtual hard particle propagator goes on-shell (referred to as the potential region) and the graviton propagator with principal value. However, it does not receive contributions from the part of the loop integrals in which the virtual graviton goes on-shell (known as the radiation-reaction region). Unfortunately, we lack a fundamental understanding of why the classical limit of the quantum soft factor does not include contributions from radiation reactions.
For a 2 → 2 scattering process with a large impact parameter or low momentum transfer, the contribution of radiation reaction to the quantum soft factor turns out to be suppressed compared to the contribution from the potential region (classical), as discussed in [1]. This has also been established in [59] using the KMOC formalism [60,61] after substituting the final momenta of the scattered particles in terms of initial momenta and the perturbatively computed momentum impulse in terms of initial scattering data and specified interaction. However, for hard scattering (small impact parameter scattering), the reason why the radiation reaction contribution to the quantum soft factor does not affect the classical waveform in the classical limit has not yet been resolved. We believe that a generalization of recent investigations into deriving classical gravitational waveforms from Eikonal exponentiation in [62][63][64][65][66][67] could potentially resolve this puzzle even for hard scattering. It would be interesting to explore whether there are any observable consequences resulting from the contribution of radiation reaction in the quantum soft factor.
There have been many applications of soft theorems both in the context of scattering amplitudes and in relation to gravitational memory. For instance, the universal characteristic of Weinberg's soft theorem imposes an infinite hierarchy of constraints on the linear momentum impulse within the KMOC formalism, as derived in [68]. It is also anticipated that the universality of the ln ω soft theorem should impose an infinite hierarchy of constraints on the angular momentum impulse within the KMOC formalism. 13 As discussed in Section 1, while the soft theorems alone cannot impose non-trivial constraints on the quantum theory in the UV, combining the results of the soft theorems with certain physical assumptions about scattering amplitudes, such as analyticity, unitarity and crossing grants them the ability to constrain the UV quantum theory. As an illustration of this concept, the article [69] derived non-perturbative bounds on the a-anomaly coefficient of the UV conformal field theory (CFT), the deformation of which leads to a massive QFT along the renormalization group (RG) flow. These bounds were established by incorporating the constraints from the double soft dilaton theorem within the framework of non-perturbative S-matrix bootstrap. For numerous astrophysical scattering events, the order of magnitude of gravitational tail memory follows from the ln ω soft theorem has been estimated in [16,20]. In certain classical scattering scenarios, the gravitational waveform resulting from the ln ω soft theorem has also been derived in [70][71][72], carrying observable consequences in the present era of gravitational wave physics.

A Intermediate steps in deriving soft photon theorem
Let us evaluate the following expression appears in (3.40) after substituting the vertex (3.12) and using the identities (2.19), (2.20) and (3.25), Now using the identities (2.11),(2.12) and (2.21) we simplify the above expression and get Let us define a specific tensor structure which will appear together in all the computations in section-3.3 , which is also the non-universal contribution to the tree-level subleading soft photon theorem as derived in [5] N i ρσ (−p i ) ≡ − i 8 With the above definition, the expression of ζ 1µ can be written in the following compact form Using the identity in (2.12) and on-shell condition (2.21) the above expression reduces to The expression below appears in (3.46) and can be evaluated analogous to the evaluation of ζ 1µ . The final result reads Let us evaluate the following expression appears in (3.32), (3.40), (3.44) and (3.46). After substituting the vertex from (3.12) and using the identities (2.19), (2.20) and (3.25) we get e j ℓ σ B j νσ (−p j )Ξ j (−p j ) + O(ℓℓ) . (A.8) Using the identities (2.11),(2.12) and the on-shell condition (2.21), the above expression reduces to ζ 4ν = e j ϵ T j (2p j + ℓ) ν + 4ℓ σ N j νσ (−p j ) + O(ℓℓ) . (A. 9) B Amputated Green's function involving single graviton Following the covariantization prescription described in section-4.1, we compute the amputated Green's function involving N number of massive spinning particles and one off-shell graviton, where the graviton is not attached to any external spinning particle leg. The resulting expression is given by Above Γ (N ) is defined after stripping out the momentum conserving delta function (2π) 4 δ (4) (p 1 + · · ·+p N ) from the expression of Γ (N ) , i.e. Γ (N ) ≡ (2π) 4 δ (4) (p 1 +· · ·+p N ) Γ (N ) . Now starting from the above covariantized expression, the goal is to express Γ (N +1) µν (ℓ) in terms of some operator operating on Γ (N ) up to linear order in ℓ. To do that we mostly follow the analysis of [3] with the only difference being h µν (ℓ) is an off-shell graviton so we can not impose traceless or transverse condition. Instead in de Donder gauge, we use ℓ µ h µν (ℓ) = 1 2 ℓ ν h(ℓ) in the intermediate stages of calculation.