§ 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 <cit.>, 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.[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 <cit.>. 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.] 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 <cit.>. 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 <cit.>. However, this issue has been successfully resolved in <cit.> by directly working in $D=4$ and carefully analyzing the possible non-analytic structures around $\omega=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\omega$, where $\omega$ 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 $\omega^0$, as the energy approaches zero ($\omega\rightarrow 0$). The existence of the $\ln\omega$ soft graviton theorem has been confirmed in <cit.> with perfect agreement with the result of <cit.> 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 <cit.>. 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 <cit.>. 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  $\lbrace p_i , \Sigma_i,\epsilon_i \rbrace $ for $i=1,2,\cdots ,N$ ($N\geq 4$) and one low-energy (soft) graviton[The graviton is the particle created by operating metric fluctuation field $h_{\mu\nu}(x)$ on vacuum, where the classical metric fluctuation is defined by $h_{\mu\nu}(x)=\f{1}{2\sqrt{8\pi G}}\left(g_{\mu\nu}(x)-\eta_{\mu\nu}\right)$.] with momentum and polarization $k,\varepsilon$, and denote this scattering amplitude by $\mathcal{A}^{(N+1)}\big(\lbrace \epsilon_i,p_i ,\Sigma_i\rbrace ,\varepsilon ,k\big)$. 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 $\omega$, so according to our convention $k^\mu=-\omega \mathbf{n}^\mu$ with $\mathbf{n}^\mu$ 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

\begin{align}
\mathcal{A}^{(N+1)}\big(\lbrace \epsilon_i,p_i ,\Sigma_i\rbrace ,\varepsilon ,k\big)=\sqrt{8\pi G}&\sum_{i=1}^{N}\epsilon_{i,\alpha} \left[(\mathbb{S}^{gr}_{\text{tree}})^{\ \alpha}_{\beta}+(\mathbb{S}^{gr}_{\text{1-loop}})^{\ \alpha}_{\beta}+(\mathbb{S}^{gr}_{\text{2-loop}})^{\ \alpha}_{\beta}+\cdots\right]\nn\\
&\times\mathcal{A}_{(i)}^{(N)\beta}\left(p_i\right)\ , \label{Soft_theorem}
\end{align}

where $\mathcal{A}_{(i)}^{(N)\beta}(p_i)$ represents the $i$'th particle polarisation $(\epsilon_{i\beta})$ stripped $N$-particle amplitude $\mathcal{A}^{(N)}\big(\lbrace \epsilon_i,p_i ,\Sigma_i\rbrace\big)$,  which is defined by the following relation
𝒜^(N)({ϵ_i,p_i ,Σ_i})≡ϵ_iβ 𝒜_(i)^(N)β(p_i) .
In (<ref>), the expression of tree level “soft factor”[The actual tree level soft factor should be think of the  expression (<ref>) with a sum over hard particles from $i=1,2,\cdots, N$. Throughout the whole paper, we follow the same terminology “soft factor” referring to the soft factor expression without hard particle sum.] for single soft graviton emission reads <cit.>
(𝕊^gr_tree)^ α_β   =    ε_μνp_i^μp_i^νp_i·kδ^α_β+ε_μνp_i^μk_ρp_i·k({p_i^νp_iρ-p_i^ρp_iν}δ^α_β+(Σ_i^ρν)_β^ α)

+𝒪(ω^n, n≥1) .
In $D=4$ analyzing tree level scattering amplitudes for effective field theory the non-universal sub-subleading soft graviton factor at order $\omega$ has also been derived in <cit.>. Reference <cit.> 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 (one-particle irreducible) vertex involving two hard particles and a soft graviton. By extending the analysis of  <cit.>, it becomes evident that a complete soft factorization is not achievable at order $\omega^n$ for $n\geq 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 <cit.>. The generalization of the tree-level soft factor $\mathbb{S}^{gr}_{\text{tree}}$ for multiple soft graviton emissions up to subleading order can be found in <cit.>.

The one-loop contribution to the “soft factor” for single soft graviton emission in (<ref>) reads
(𝕊^gr_1-loop)^ α_β   =     K_phase^reg ε_μνp_i^μp_i^νp_i·kδ^α_β+ε_μνp_i^μk_ρp_i·k {p_i^νK_gr^regp_iρ-p_i^ρK_gr^regp_iν}δ^α_β

+ K_phase^reg ε_μνp_i^μk_ρp_i·k({p_i^νp_iρ-p_i^ρp_iν}δ^α_β+(Σ_i^ρν)_β^ α)

+ε_μν k_ρk_σp_i·k{p_i^μK_gr^regp_iρ-p_i^ρK_gr^regp_iμ}({p_i^νp_iσ-p_i^σp_iν}δ^α_β+(Σ_i^σν)_β^ α)

+ 𝒪(ω^n,n≥0)  + 𝒪(ω^nlnω,n≥2) ,

\begin{align}
K^{reg}_{gr}\ =\ \f{i}{2}\ (8\pi G)\sum_{\ell=1}^{N}\  \sum_{\substack{j=1\\j\neq \ell}}^{N}\Big{\lbrace} (p_{\ell}.p_{j})^{2}-\f{1}{2}p_{\ell}^{2}p_{j}^{2}\Big{\rbrace}\ \int_{\omega}^{\Lambda} \f{d^{4}\ell}{(2\pi)^{4}} \ \f{1}{\ell^{2}-i\epsilon}\ \f{1}{(p_{\ell}\cdot \ell +i\epsilon)\ (p_{j}\cdot \ell - i\epsilon)}\nn\\
\simeq  -iG\ (\ln\omega)\sum_{\ell=1}^{N} \sum_{\substack{j=1\\j\neq \ell}}^{N} \ \f{\Big{\lbrace} (p_{\ell}.p_{j})^{2}-\f{1}{2}p_{\ell}^{2}p_{j}^{2}\Big{\rbrace}}{\sqrt{(p_{\ell}.p_{j})^{2}-p_{\ell}^{2}p_{j}^{2}}} \Bigg{\lbrace}\delta_{\eta_{\ell}\eta_{j},1}
-\f{i}{2\pi}\ln\Bigg(\f{p_{\ell}.p_{j}+\sqrt{(p_{\ell}.p_{j})^{2}-p_{\ell}^{2}p_{j}^{2}}}{p_{\ell}.p_{j}-\sqrt{(p_{\ell}.p_{j})^{2}-p_{\ell}^{2}p_{j}^{2}}}\Bigg)\Bigg{\rbrace}\ , \label{K_gr_reg}
\end{align}

K_phase^reg    =    i (8πG) ∑_j=1^N(p_j.k)^2∫_R^-1^ωd^4ℓ(2π)^4 1ℓ^2-iϵ1k.ℓ+iϵ1p_j.ℓ-iϵ

≃    -2iG (lnω)  [ ∑_j=1
η_j=-1^Np_j.k - i2π∑_j=1^Np_j.k ln(p_j^2(p_j.𝐧)^2)] .
In the expressions (<ref>) and (<ref>), $\eta_j=+1$ if $j$'th particle is ingoing and $\eta_j=-1$ if $j$'th particle is outgoing. Under the sign $\simeq $ we only keep the logarithmic contributions after performing the integrations following <cit.>. The upper limit of the integration $\Lambda$ in (<ref>) represents the order of the energy of hard particles and  the lower limit of the integration $R^{-1}$ in (<ref>) represents the order of the energy resolution of the detector. The first line of (<ref>) contains the $\mathcal{O}(\ln\omega)$ soft factor which has been derived in  <cit.> as an one-loop exact result, analyzing one-loop gravitational S-matrices in the theory of scalar coupled to gravity. In <cit.> a correction to $\mathcal{O}(\ln\omega)$ 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 $\mathcal{O}(\ln\omega)$ 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 $\mathcal{O}(\ln\omega)$ soft factor, while also extending the infrared divergence factorization prescription proposed in <cit.> to encompass a broad range of quantum gravity theories. In this paper we also conjecture the order $\omega\ln\omega$ soft factor given in the second and third lines of (<ref>) 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 $\mathcal{O}(\omega\ln\omega)$ soft factor has already been derived in <cit.> in the name of spin-dependent classical soft graviton theorem which provides evidence on the correctness of the conjecture. Note that the tree level subleading soft theorem result at order $\omega^0$ in (<ref>) 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) <cit.>.

The two-loop contribution to the “soft factor” for single soft graviton emission in (<ref>) reads
(𝕊^gr_2-loop)^ α_β   =      12{K_phase^reg}^2 ε_μνp_i^μp_i^νp_i·kδ^α_β+ K_phase^regε_μνp_i^μk_ρp_i·k(p_i^νK_gr^regp_iρ-p_i^ρK_gr^regp_iν)δ^α_β

+ 12∑_i=1^Nε_μνk_ρk_σp_i·k (p_i^μK_gr^regp_iρ-p_i^ρK_gr^regp_iμ)(p_i^νK_gr^regp_iσ-p_i^σK_gr^regp_iν) δ^α_β

+ 𝒪(ω^n(lnω)^2, n≥2) + 𝒪(ω^nlnω, n≥1) +𝒪(ω^n, n≥1) .
The above result was conjectured in <cit.> 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 $\mathcal{O}\big(\omega(\ln\omega)^2\big)$, was derived under the name of the sub-subleading classical soft graviton theorem in <cit.>. This classical derivation offers substantial evidence supporting the validity of the above two-loop soft factor. Note that the $\mathcal{O}(\omega \ln \omega)$ soft factor in the second and third lines of (<ref>) 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 $\omega$ approaches zero, behaves like $\omega^{n-1}(\ln\omega)^n$ relating it to tree level $N$-particle amplitude. The general structure of the order $\omega^{n-1}(\ln\omega)^n$ soft graviton theorem is provided in <cit.>.

Note that the “tree”, “1-loop” , ... subscripts in the soft factors in the expression (<ref>) 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 $\mathbb{S}^{gr}_{\text{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\geq 1$. Similarly $\mathbb{S}^{gr}_{\text{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\geq 2$, and $\mathbb{S}^{gr}_{\text{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\geq 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 $\lbrace e_i \rbrace$ and study scattering amplitude involving one soft photon emission with polarization and momentum $(\varepsilon , k )$. So in this case the soft expansion of $(N+1)$ particle amplitude takes the following form

\begin{align}
\mathcal{A}^{(N+1)}\big(\lbrace \epsilon_i,p_i ,e_i,\Sigma_i\rbrace ,\varepsilon ,k\big)=&\sum_{i=1}^{N}\epsilon_{i,\alpha} \left[(\mathbb{S}^{em}_{\text{tree}})^{\ \alpha}_{\beta}+(\mathbb{S}^{em}_{\text{1-loop}})^{\ \alpha}_{\beta}+(\mathbb{S}^{em}_{\text{2-loop}})^{\ \alpha}_{\beta}+\cdots\right]\nn\\
&\times\mathcal{A}_{(i)}^{(N)\beta}\left(p_i\right)\ , \label{soft_photon_theorem}
\end{align}

where $\mathcal{A}_{(i)}^{(N)\beta}\left(p_i\right)$ is defined through the relation (<ref>). In (<ref>), the expression of tree level “soft factor” for single soft photon emission reads <cit.>
(𝕊^em_tree)^ α_β=      e_iε_μp_i^μp_i·kδ^α_β+e_iε_μ k_ρp_i·k({p_i^μp_iρ-p_i^ρp_iμ}δ^α_β+(𝒩_i^ρμ(-p_i))_β^ α)

where the generic expression for the non-universal term $ \mathcal{N}_{i}^{\rho\sigma}(-p_i)$ contributing to the subleading soft photon theorem has been derived in <cit.> and its explicit form in provided in (<ref>). $ \mathcal{N}_{i}^{\rho\mu}(-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  <cit.> it can be argued that a complete soft factorization of order $\omega^n$ for $n\geq 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 <cit.>. The generalization of the tree-level soft factor $\mathbb{S}^{em}_{\text{tree}}$ for multiple soft photon emissions up to subleading order can be found in <cit.>.

The one-loop contribution to the “soft factor” for single soft photon emission in (<ref>) reads
(𝕊^em_1-loop)^ α_β   =    e_iε_μ k_ρ p_i·k {p_i^μK_em^regp_iρ-p_i^ρK_em^regp_iμ}δ^α_β

+𝒪(ω^nlnω, n≥1) +𝒪(ω^n, n≥0) ,
K_em^reg= i2∑_ℓ=1^N  ∑_j=1
j≠ℓ^Ne_ℓe_j(p_ℓ·p_j) ∫_ω^Λd^4ℓ(2π)^4  1ℓ^2-iϵ 1(p_ℓ·ℓ+iϵ) (p_j·ℓ- iϵ)

≃    -i2∑_ℓ=1^N  ∑_j=1
j≠ℓ^Ne_ℓe_j4π (lnω)  p_ℓ·p_j√((p_ℓ.p_j)^2-p_ℓ^2p_j^2) {δ_η_ℓη_j,1-i2πln(p_ℓ.p_j+√((p_ℓ.p_j)^2-p_ℓ^2p_j^2)p_ℓ.p_j-√((p_ℓ.p_j)^2-p_ℓ^2p_j^2))} .

The $\mathcal{O}(\ln\omega)$ soft factor for single soft photon emission in (<ref>) has been derived in  <cit.> as an one-loop exact result, analyzing one-loop S-matrices in the theory of minimally coupled scalar QED. There a correction to $\mathcal{O}(\ln\omega)$ soft photon factor due to gravitational interaction has also been derived. In this article we re-derive the $\mathcal{O}(\ln\omega)$ 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 $\mathcal{O}(\ln\omega)$ soft factor, while also extending the infrared divergence factorization prescription proposed in <cit.> to encompass a broad range effective field theories for charged objects. The generalization of the one-loop soft factor $\mathbb{S}^{em}_{\text{1-loop}}$ for multiple soft photon emissions up to subleading order can be found in the section-(3.5) of <cit.>, and the final result has been provided in (<ref>).

The two-loop contribution to the “soft factor” for single soft photon emission in (<ref>) reads
(𝕊^em_2-loop)^ α_β   =    e_i2ε_μk_ρp_i·k (p_i^μK_em^regp_iρ - p_i^ρK_em^regp_iμ ) (  k_σK_em^regp_iσ)δ^α_β

+ (lnω)^2 e_i( ε^μ k^ρ-ε^ρk^μ) p_iμ𝒞^reg_ρ(q_i,p_i;{e_j},{p_j})δ^α_β

+ 𝒪(ω^n(lnω)^2, n≥2)+ 𝒪(ω^nlnω, n≥1) +𝒪(ω^n, n≥1) ,

=   -∑_j=1
ℓ≠i^Ne_i^2e_je_ℓ4 {p_i.p_jδ_ρ^κ-p_i^κp_jρ}p_iσ{ℐ(p_i,p_ℓ)×p_i.p_ℓ} ^2 ℐ(p_i,p_j)p_i^σp_i^κ

ℓ≠j^Ne_ie_j^2e_ℓ4 {p_i.p_jδ_ρ^κ-p_i^κp_jρ}p_jσ{ℐ(p_j,p_ℓ)×p_j.p_ℓ}^2 ℐ(p_i,p_j)p_i^σp_i^κ ,
ℐ(p_i,p_j)= -14π 1√((p_i.p_j)^2-p_i^2p_j^2){δ_η_iη_j,1-i2πln(p_i.p_j+√((p_i.p_j)^2-p_i^2p_j^2)p_i.p_j-√((p_i.p_j)^2-p_i^2p_j^2))}.
The $\mathcal{O}\left(\omega(\ln\omega)^2\right)$ soft factor for single soft photon emission in (<ref>) has been derived in the section-4 of <cit.> 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 $\omega$ approaches zero, behaves like $\omega^{n-1}(\ln\omega)^n$ and it relates to the tree level $N$-particle amplitude. The general structure of the order $\omega^{n-1}(\ln\omega)^n$ soft photon theorem is provided in <cit.>.

The rest of the paper is organized as follows:
In section-<ref>, 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-<ref>, 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-<ref>, 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-<ref> and <ref>, we discuss the possible generalizations of our derivations to higher orders.
In section-<ref>, 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, \ldots$ as Lorentz indices for the tangent space, ranging from 0 to 3. The Latin alphabets starting from $i, j, k, \ell, \ldots$ are employed as indices for identifying individual hard particles, ranging from 1 to N. The first few Greek alphabets $\alpha, \beta, \gamma, \delta, \ldots$ are used as polarization indices for spinning particles on the tangent space, while the Greek alphabets beginning with $\lambda, \mu, \nu, \rho, \sigma, \tau, \ldots$ serve as curved space indices, ranging from 0 to 3. In section-<ref>, where we derive the soft photon theorem solely under electromagnetic interaction in a flat background, we employ both $a, b, c, d, \ldots$ and  $\lambda, \mu, \nu, \rho, \sigma, \tau, \ldots$ as flat space Lorentz indices.

Metric and unit conventions: In our convention four dimensional Minkowski metric is $\eta_{ab}=\text{diag}(-1,+1,+1,+1)$. We work in the unit where speed of light $c=1$ and Planck constant $\hbar=1$ but keep the gravitational constant $G$ explicit. We define $\kappa\equiv \sqrt{8\pi 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  $\lbrace e_i, p_i , \Sigma_i,\epsilon_i \rbrace $ for $i=1,2,\cdots ,N$ and one low-energy (soft) outgoing photon/graviton with momenta and polarization $k,\varepsilon$, and denote this scattering amplitude by $\mathcal{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 $\omega$ so that $k^{\mu}=-\omega \mathbf{n}^{\mu}$ where $\mathbf{n}^\mu$ being the null vector whose spatial part denotes the direction of soft photon/graviton emission. Here we are only interested to evaluate $\mathcal{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 ($\omega<<|p_i|$) of $\mathcal{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  $\mathcal{A}^{(N)}$.[Note that the soft limit can also be defined covariantly by demanding $\Big{|}\f{p_i.k}{p_i.p_j}\Big{|}<<1$ for all $i,j=1,\cdots ,N$.] Note that both the scattering amplitudes $\mathcal{A}^{(N)}$ and $\mathcal{A}^{(N+1)}$ are distributions in momenta as  $\mathcal{A}^{(N)}$ contains momentum conserving delta function $\delta^{(4)}\big(p_1+p_2+\cdots +p_N\big)$ and $\mathcal{A}^{(N+1)}$ contains momentum conserving delta function $\delta^{(4)}\big(p_1+p_2+\cdots +p_N+k\big)$. 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 $\mathcal{A}^{(N+1)}$ and $\mathcal{A}^{(N)}$ after full/partial cancellation of the IR divergent contributions as we explain in later sections.

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 <cit.>.

Handling IR divergences in the derivation of soft theorem:
The traditional S-matrix 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 <cit.>, 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 <cit.>. A generalization of Grammer-Yennie prescription for perturbative QCD and quantum gravity can be found in <cit.> and  <cit.> respectively.

In the derivation of soft photon theorem, Grammer-Yennie prescription helps to factor out IR divergences from both the amplitudes $\mathcal{A}^{(N+1)}$ and $\mathcal{A}^{(N)}$ in the following way
𝒜^(N)=exp{K_em} 𝒜^(N)_IR-finite , 𝒜^(N+1)=exp{K_em} 𝒜^(N+1)_IR-finite .
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 (<ref>). Basically Grammer-Yennie prescription provides a systematic procedure to compute IR-finite parts perturbatively for both the amplitudes. When the soft factor $\mathbb{S}^{em}$ is a multiplicative function instead of differential operator, we get
𝒜^(N+1)=𝕊^em  𝒜^(N) ⟹ 𝒜^(N+1)_IR-finite=𝕊^em  𝒜^(N)_IR-finite .
Hence to derive the $\mathcal{O}(\ln\omega)$ and $\mathcal{O}\left(\omega(\ln\omega)^2\right)$ soft factors in (<ref>) and (<ref>) we can directly start from one and two-loop contributions of IR finite amplitude $\mathcal{A}^{(N+1)}_{\text{IR-finite}}$ and perform soft expansion.

In the derivation of soft graviton theorem, $\mathcal{A}^{(N+1)}$ contains some extra divergent factors relative to $\mathcal{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 <cit.> would be
𝒜^(N)=exp{K_gr} 𝒜^(N)_IR-finite , 𝒜^(N+1)=exp{K_gr+K_phase} 𝒜^(N+1)_IR-finite ,
where the IR-divergent expressions of $K_{gr}$ and $K_{phase}$ are given in (<ref>) and (<ref>). 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 $\mathbb{S}^{gr}$ is a multiplicative function instead of differential operator, we get
𝒜^(N+1)=𝕊^gr  𝒜^(N) ⟹ 𝒜^(N+1)_IR-finite=exp{-K_phase} 𝕊^gr  𝒜^(N)_IR-finite .
Hence to derive an unambiguous soft factor by analyzing $\mathcal{A}^{(N+1)}_{\text{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 $\mathcal{O}(\ln\omega)$ and $\mathcal{O}\left(\omega(\ln\omega)^2\right)$ soft factors in (<ref>) and (<ref>) respectively by analyzing one and two loop IR finite amplitudes. But if we want to derive the $\mathcal{O}(\omega\ln\omega)$ soft factor in the second and third lines of (<ref>) we need to deal with the following additional subtleties:

* Since the order $\omega\ln\omega$ soft factor in (<ref>) is a differential operator we can not really commute the soft factor and infrared divergent exponential to get the second equation of (<ref>). Hence to derive the $\mathcal{O}(\omega\ln\omega)$ soft factor we have to start with the full divergent scattering amplitude $\mathcal{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.

* Note that the momentum conserving delta function associated with $\mathcal{A}^{(N+1)}$ is
$\delta^{(4)} \left(\sum\limits_{i=1}^N p_i+k\right)$. On the other hand the momentum conserving delta function associated with $\mathcal{A}^{(N)}$ is $\delta^{(4)} \left(\sum\limits_{i=1}^N p_i\right)$. Now Taylor series expansion of the first delta function around small $\omega$ produces a correction of order $\omega$. This correction, when multiplied with the $\mathcal{O}(\ln\omega)$ soft factor, yields an additional factor of order $\omega\ln\omega$ 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  $\omega\ln\omega$ soft graviton factor in (<ref>) for future study and focus on deriving the order $\ln\omega$ 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 one-particle 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 $\Phi_{\alpha}(x)$ denotes the set of all massive fields in real representation, present in the 1PI effective action[In $\Phi_{\alpha}(x)$ we also can include massless finite energy particles e.g. hard photon or graviton, in that case our 1PI effective action should be thought of as tree level action.] which transforms in a reducible representation of Lorentz group $SO(1,3)$ in the following way,
SO(1,3):        x^a →x^'a=Λ^a _b  x^b = (δ^a_b +λ^a _b)x^b +𝒪(λ^2)

Φ_α(x)→Φ'_α(x)= [exp{-i2λ_abΣ^ab}]_α^ β Φ_β(Λ^-1x) ,
where $\lambda_{ab}=-\lambda_{ba}$ is the infinitesimal Lorentz transformation parameter and $\Sigma^{ab}$ is the spin angular momentum generator of $SO(1,3)$ transformation in the real reducible representation. The subscript index $\alpha$ 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 $\Phi(x)$ transforms in the following way,
U(1)_EM: Φ_α(x)→Φ'_α(x)= [exp{i 𝒬 θ}]_α^ β Φ_β(x) ,
where $\theta$ is the parameter of global $U(1)_{EM}$ transformation and $\mathcal{Q}$ is the generator of global $U(1)_{EM}$ transformation in the real representation of $\Phi(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 <cit.>, 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 $\mathcal{Q}$ is the generator of $SO(2)$ transformation. In the set of fields denoted by $\Phi_{\alpha}(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 $\mathcal{Q}$ will be zero.

Let us start with the general form of the quadratic part of the massive particle 1PI effective action[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 $\mathcal{K}(q)$ may not be polynomially expandable around $q^\mu =0$, which is the key assumption for the validity of the covariantization prescription discussed below. ]
S^(2)   =   12∫d^4 q_1(2π)^4d^4 q_2(2π)^4 (2π)^4δ^(4)(q_1+q_2) Φ_α(q_1)𝒦^αβ(q_2)Φ_β(q_2)

=   12∫d^4 q_1(2π)^4d^4 q_2(2π)^4 (2π)^4δ^(4)(q_1+q_2) Φ^T(q_1)𝒦(q_2)Φ(q_2) ,
where $\mathcal{K}(q)$ is the renormalized momentum space kinetic operator which satisfy the following condition:
𝒦^αβ(q)=± 𝒦^βα(-q) ,

𝒦(q)=± [𝒦(-q)]^T .
In the second lines of (<ref>) and (<ref>) 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 $\Phi(x)$ 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 (<ref>) becomes
Δ^i_αβ(q) = i[𝒦_i^-1(q)]_αβ≡(q^2+m_i^2-iϵ)^-1 Ξ_iαβ(q) ,
where $\mathcal{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 $\Xi_i(q)$ as the residue of the pole of the propagator for $i$-th particle. The relation between $\mathcal{K}_i$ and $\Xi_i$ and their momentum derivatives satisfy the following relations in the index free notation, which will be useful for later computation <cit.>
𝒦_i(q)Ξ_i(q)=i(q^2+m_i^2 -iϵ) ,

𝒦_i(q)q^aΞ_i(q)=-𝒦_i(q)Ξ_i(q)q^a+2iq_a ,

^2 𝒦_i(q)q^a q^bΞ_i(q)=-𝒦_i(q)q^aΞ_i(q)q^b-𝒦_i(q)q^bΞ_i(q)q^a-𝒦_i(q)^2 Ξ_i(q)q^aq^b+2iη_ab ,

Ξ_i(q)𝒦_i(q)=i(q^2+m_i^2-iϵ) ,

Ξ_i(q)q^a𝒦_i(q)=-Ξ_i(q)𝒦_i(q)q^a+2iq_a ,

^2 Ξ_i(q)q^a q^b𝒦_i(q)=-Ξ_i(q)q^a𝒦_i(q)q^b-Ξ_i(q)q^b𝒦_i(q)q^a-Ξ_i(q)^2 𝒦_i(q)q^aq^b+2iη_ab .
The Lorentz covariance of $\mathcal{K}_i$ and $\Xi_i$ implies the following two relations
(Σ_i^ab)^T𝒦_i(q)=-𝒦_i(q)Σ_i^ab+q^a 𝒦_i(q)q_b-q^b 𝒦_i(q)q_a ,

Σ_i^ab Ξ_i(q)=-Ξ_i(q) (Σ_i^ab)^T-q^a Ξ_i(q)q_b+q^b Ξ_i(q)q_a ,
where $\Sigma_i$ is the spin angular momentum generator for $i$-th component field inside $\Phi(x)$. Taking derivatives with respect to momenta the above expressions become
(Σ_i^ab)^T𝒦_i(q)q_c=-𝒦_i(q)q_cΣ_i^ab+q^a ^2 𝒦_i(q)q_b q_c-q^b ^2 𝒦_i(q)q_a q_c+η^ac𝒦_i(q)q_b-η^bc𝒦_i(q)q_a,

Σ_i^ab Ξ_i(q)q_c=-Ξ_i(q)q_c (Σ_i^ab)^T-q^a ^2 Ξ_i(q)q_b q_c+q^b ^2 Ξ_i(q)q_aq_c-η^acΞ_i(q)q_b+η^bcΞ_i(q)q_a .

Invariance of (<ref>) under global $U(1)_{EM}$ transformation implies
𝒬_γ ^α𝒦^γβ +𝒦^αγ𝒬_γ ^β =0 ⇒ 𝒬^T𝒦+𝒦𝒬=0 .
This also imposes constraint on the numerator of the propagator which reads
𝒬Ξ+Ξ𝒬^T=0 .
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 $\epsilon_i(q_i)$ it satisfies
𝒦_i^αβ(q_i)ϵ_iβ(q_i)=0  ⇒ 𝒦_i(q_i)ϵ_i(q_i)=0 and ϵ_i^T(q_i)𝒦_i^T(q_i)=0 .

§ 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 <cit.>, the order $\ln\omega$ soft factor has been derived analyzing one-loop 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\omega$ soft factor derived in  <cit.> 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 <cit.>, the quadratic action $S^{(2)}$ in (<ref>) 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
Δ^F_μν(ℓ)=-iℓ^2-iϵ η_μν .
In position space the kinetic operator in (<ref>) contains derivatives over the field $\Phi_{\beta}(x)$, which have to replace by covariant derivatives under covariantization in presence of photon field $A_{\mu}(x)$. For example in the case of one and two derivatives the covariantization rule in position space becomes
_μΦ_β→ D_μΦ_β      = (_μδ_β^γ-i𝒬_β ^γA_μ)Φ_γ ,

_μ_νΦ_β→ D_(μ D_ν)Φ_β      = [_(μ δ_β^γ-i𝒬_β ^γA_(μ][_ν) δ_γ^δ-i𝒬_γ ^δA_ν)]Φ_δ

=_μ_νΦ_β-i𝒬_β ^γ(A_μ_ν+A_ν_μ)Φ_γ

-i2𝒬_β ^γ(_μA_ν+_νA_μ)Φ_γ-A_μA_ν𝒬_β ^γ 𝒬_γ ^δΦ_δ .
Above we use the following symmetrization convention $E_{(\mu}F_{\nu)}=\f{1}{2}\big(E_\mu F_\nu +E_\nu F_\mu\big)$ for two vectors $E$ and $F$. For determining minimal interaction vertices $\Gamma^{(3)}$ involving two massive spinning particles and one photon up to one derivative on photon field, and $\Gamma^{(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 (<ref>):
S^(3)   =   12∫d^4 q_1(2π)^4d^4 q_2(2π)^4d^4ℓ(2π)^4 (2π)^4δ^(4)(q_1+q_2+ℓ)Φ_α(q_1)[-A^μ(ℓ)𝒦^αγ(q_2)q_2^μ𝒬_γ ^β

-14(ℓ^μA^ν(ℓ)+ℓ^νA^μ(ℓ))^2 𝒦^αγ(q_2)q_2^μq_2^ν𝒬_γ ^β + 𝒪(ℓ^2)]Φ_β(q_2) ,
S^(4)   =   12∫d^4 q_1(2π)^4d^4 q_2(2π)^4d^4ℓ_1(2π)^4d^4ℓ_2(2π)^4 (2π)^4δ^(4)(q_1+q_2+ℓ_1+ℓ_2)

Φ_α(q_1)[ 12 A^μ(ℓ_1) A^ν(ℓ_2)^2 𝒦^αδ(q_2)q_2^μq_2^ν𝒬_δ ^γ𝒬_γ ^β + 𝒪(ℓ_1,ℓ_2)]Φ_β(q_2) .
Above $A_\mu(\ell)$ is the Fourier transform of gauge field defined through the following relation: $A_\mu(x)\equiv \int \f{d^4\ell}{(2\pi)^4}\ e^{i\ell\cdot x}A_\mu(\ell)$. 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
S̅^(3)= 12       ∫d^4q_1(2π)^4d^4q_2(2π)^4d^4ℓ(2π)^4 (2π)^4δ^(4)(q_1+q_2+ℓ)

Φ_α(q_1) [F_μν(ℓ)  ℬ^αβ,μν(q_2) +𝒪(ℓ^2)] Φ_β(q_2) .
F_μν(ℓ)=i[ℓ_μA_ν(ℓ)-ℓ_νA_μ(ℓ)] ,
and $\mathcal{B}(q_2)$ satisfies the following relations
𝒬_γ ^αℬ^γβ,μν(q_2)+ℬ^αγ,μν(q_2)𝒬_γ ^β=0       ⇒ 𝒬^Tℬ^μν(q_2)+ℬ^μν(q_2)𝒬=0 ,

ℬ^αβ,μν(q_2) = -ℬ^αβ,νμ(q_2)       ⇒ ℬ^μν(q_2)=-ℬ^νμ(q_2) ,

ℬ^αβ,μν(q_2) = ± ℬ^βα,μν(-q_1-ℓ)       ⇒ℬ^μν(q_2) = ± (ℬ^μν(-q_1-ℓ))^T  .
In the last equation above $+$ sign is for Grassmannian even field and $-$ sign is for Grassmannian odd field. Again during the derivation we consider $\Phi$ 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 (<ref>) and (<ref>), the interaction vertex describing two ingoing spinning particles with momenta $q$ and $-(q+\ell)$, polarization index $\alpha$ and $\beta$, and one ingoing photon with momentum $\ell$ becomes

=    i2[ 𝒦^αγ(-q-ℓ)q^μ𝒬_γ ^β-12ℓ^ν^2𝒦^αγ(-q-ℓ)q^μq^ν𝒬_γ ^β -2iℓ^νℬ^αβ, _μν(-q-ℓ)

-𝒦^βγ(q)q^μ𝒬_γ ^α-12ℓ^ν^2𝒦^βγ(q)q^μq^ν𝒬_γ ^α -2iℓ^νℬ^βα, _μν(q) +𝒪(ℓ^2)] ,
Expanding in small $\ell$ limit and using (<ref>),(<ref>), (<ref>) the above vertex reduce to the following polarization index suppressed form
=    i[ 𝒦(-q)q^μ𝒬+12ℓ^ν^2 𝒦(-q) q^μq^ν𝒬-2iℓ^νℬ_μν(-q) + 𝒪(ℓ^2)].
Similarly starting from (<ref>) the polarization index suppressed four point interaction vertex involving two incoming spinning particles with momenta $q$ and $-(q+\ell_1+\ell_2)$ and two incoming photons with momenta $\ell_1$ and $\ell_2$ becomes
Γ^(4)_μν(q,-q-ℓ_1-ℓ_2,ℓ_1, ℓ_2)   =    i^2 𝒦(-q)q^μq^ν𝒬𝒬+ 𝒪(ℓ_1,ℓ_2) .
We denote the scattering amplitude describing $N$-number of spinning hard particle scattering in massive EFT by $\Gamma^{(N)}$, which can be expressed as a polarization tensor contracted form in the following way
Γ^(N)    =    {∏_i=1^Nϵ_α_i(p_i) } Γ^(N)α_1α_2⋯α_N

=     ϵ_α_i(p_i)ϵ_α_j(p_j) Γ_(ij)^(N)α_iα_j(p_i,p_j)

=     ϵ_α_i(p_i)Γ_(i)^(N)α_i(p_i) ,
Γ_(ij)^(N)α_iα_j(p_i,p_j)   ≡   {∏_k=1
k≠i,j^Nϵ_α_k(p_k) } Γ^(N)α_1α_2⋯α_N ,

Γ_(i)^(N)α_i(p_i)    ≡    {∏_j=1
j≠i^Nϵ_α_j(p_j) } Γ^(N)α_1α_2⋯α_N.
Above $\epsilon_{\alpha_i}(p_i)$ denotes the polarization tensor for $i$-th spinning particle with momentum $p_i$. Here we also should remember that $\Gamma^{(N)}$ is a distribution as it contains a momentum conserving delta function $\delta^{(4)}\big(p_1+p_2+\cdots +p_N\big)$. We denote a part of amputated Green's function involving $N$ number of spinning hard particles and one  photon with momentum $\ell$ by $\widetilde{\Gamma}_\mu^{(N+1)}(\ell)$, which describes sum of the contributions of the Feynman diagrams where the photon is not connected to any external leg. Since $\widetilde{\Gamma}_\mu^{(N+1)}(\ell)$ does not contain any loop involving massless particles, we can write down the following relation between $\widetilde{\Gamma}_\mu^{(N+1)}(\ell)$ and $\Gamma^{(N)}$ using the same covariantization prescription described above
Γ_μ^(N+1)α_1⋯α_N(ℓ)=-∑_i=1^N 𝒬_β_i ^α_i p_i^μ Γ^(N)α_1⋯α_i-1β_iα_i+1⋯α_N + 𝒪(ℓ)  .
In the above expression $\Gamma^{(N)}$ contains momentum conserving delta function $\delta^{(4)}\big(p_1+p_2+\cdots +p_N\big)$ and $\widetilde{\Gamma}^{(N+1)}$ contains momentum conserving delta function $\delta^{(4)}\big(p_1+p_2+\cdots +p_N+\ell\big)$.

§.§ Grammer-Yennie decomposition and IR-finite amplitudes

In this section we discuss the Grammer-Yennie decomposition introduces in <cit.> and show how it helps to factorize the IR-divergent exponential from the IR-finite part of the amplitudes as proposed in (<ref>).
In Feynman gauge we decompose  the photon propagator with momentum $\ell$ flowing from the leg $i$ to the leg $j$ for $i\neq j$
Δ_F^μν(ℓ)=-i η^μν ℓ^2 - iϵ    =    -i ℓ^2-iϵ
{K_(ij)^μν + G_(ij)^μν} ,

K_(ij)^μν = ℓ^μℓ^ν(2p_i-ℓ)·(2p_j+ℓ) (2p_i·ℓ-ℓ^2+iϵ)
(2p_j·ℓ+ℓ^2-iϵ)  ,
G_(ij)^μν = η^μν - K_(ij)^μν  .
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 (<ref>). 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 (<ref>) 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 $\ell^\mu\ell^\nu$ 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.<ref>. The LHS of Fig.<ref> after contracting with $\ell^\mu$ representing K-photon can be expressed as

This figure is a Feynman diagrammatic representation of the expression in (<ref>). Solid lines represent the massive spinning particles, dashed lines represent the ingoing virtual photon with momentum $\ell$ 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 $\mathcal{Q}_i$ and $\mathcal{Q}_i^T$ as written next to the vertices.

1q_i^2+m_i^2-iϵ Ξ_i(-q_i) ℓ^μΓ^(3)_μ(q_i,-q_i-ℓ,ℓ) Ξ_i(-q_i-ℓ) 1(q_i+ℓ)^2+m_i^2 -iϵ

=    i1q_i^2+m_i^2-iϵ1(q_i+ℓ)^2+m_i^2 -iϵ[-i(2q_i.ℓ+ℓ^2)Ξ_i(-q_i)𝒬_i^T

+Ξ_i(-q_i)𝒦_i(-q_i){ℓ^μΞ_i(-q_i)q_i^μ+12ℓ^μℓ^ν^2 Ξ_i(-q_i)q_i^μq_i^ν + 𝒪(ℓ^3)}𝒬^T_i] .
To write down the Feynman rule in the first line above for the diagram in Fig.<ref>, we follow the convention that our time arrow runs from right to left. We shall follow the same time arrow convention for all the Feynman diagrams in this article. To evaluate the first line, with the expression in (<ref>) we use the properties (<ref>),(<ref>),(<ref>) and (<ref>). Now replacing  $(2q_i.\ell +\ell^2)=\lbrace (q_i+\ell)^2+m_i^2\rbrace -\lbrace q_i^2+m_i^2\rbrace$ within the square bracket above and using the property (<ref>), the RHS of the above expression translates to

-1(q_i+ℓ)^2+m_i^2 -iϵ
{ℓ^μΞ_i(-q_i)q_i^μ+12ℓ^μℓ^ν^2 Ξ_i(-q_i)q_i^μq_i^ν + 𝒪(ℓ^3)}𝒬^T_i .
If we un-do the small $\ell$ expansion in the second line of the above expression[Instead of (<ref>), if we use the unexpanded expression (<ref>) for $\Gamma^{(3)}$ vertex, we do not need to un-do the  small $\ell$ expansion to derive the result below. Actually the relation in (<ref>) is an exact relation valid up to all order in $\ell$ expansion with any arbitrary non-minimal coupling contributing to $\Gamma^{(3)}$, as it is a direct consequence of the Ward Identity.]  and use the relation (<ref>) the  Ward identity turns out to be the following expression which has been diagrammatically represented in Fig.<ref>.
1q_i^2+m_i^2-iϵ Ξ_i(-q_i) ℓ^μΓ^(3)_μ(q_i,-q_i-ℓ,ℓ) Ξ_i(-q_i-ℓ) 1(q_i+ℓ)^2+m_i^2 -iϵ

=   𝒬_iΞ_i(-q_i-ℓ)(q_i+ℓ)^2+m_i^2-iϵ + Ξ_i(-q_i)q_i^2+m_i^2-iϵ𝒬_i^T .
Important to note that the solid blob vertices in Fig.<ref> 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 $\nu$. 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.<ref>. Using the Ward identity described in Fig.<ref> for first and third diagrams in the first line of Fig.<ref>, we find the diagrams drawn after the equality in Fig.<ref>.  Now if we can show that the sum of the contribution of the three diagrams in the last line of Fig.<ref> vanishes, then the Ward identity of Fig.<ref> is also valid in presence of an external photon line. The three diagrams in the last line of Fig.<ref> has been drawn again in Fig.<ref> and the contribution becomes
1q_i^2+m_i^2-iϵ Ξ_i(-q_i)[ 𝒬_i^T Γ^(3)_ν(q_i+ℓ,-q_i-ℓ-k, k) + ℓ^μ Γ^(4)_μν(q_i,-q_i-ℓ-k,ℓ, k)

+ Γ^(3)_ν(q_i, -q_i-k,k) 𝒬_i ] Ξ_i(-q_i-ℓ-k)1(q_i+ℓ+k)^2 +m_i^2 -iϵ .

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.<ref> we get the diagrams in second and third lines.

Substituting the expressions for the vertices from (<ref>) and (<ref>), the expression inside the square bracket of (<ref>)  turns out to be
i𝒬_i^T[ 𝒦_i(-q_i-ℓ)q_i^ν𝒬_i +12k^μ^2 𝒦_i(-q_i-ℓ)q_i^μq_i^ν𝒬_i -2ik^μℬ^i_νμ(-q_i-ℓ)]

+ iℓ^μ ^2 𝒦_i(-q_i)q_i^μq_i^ν 𝒬_i𝒬_i

+i[ 𝒦_i(-q_i)q_i^ν𝒬_i +12k^μ^2 𝒦_i(-q_i)q_i^μq_i^ν𝒬_i -2ik^μℬ^i_νμ(-q_i)]𝒬_i .

Identity involving last three Feynman diagrams in Fig.<ref> whose mathematical expression has been written in (<ref>).

Generalization of the identity in Fig.<ref> in presence of $n$-number of photons.

By Taylor expanding the first line above for small $\ell$ and only keeping terms up to linear order in $\ell$ or $k$, we can use the identities in (<ref>) and (<ref>) to show that the sum of the contribution in the three lines above vanishes up to linear order in $\ell$ or $k$. This proves the diagrammatic identity in Fig.<ref>. By utilizing the expressions of vertices from (<ref>) and (<ref>), which are given up to linear and zeroth order in $\ell$ or $k$ respectively, it may initially appear that the validity of the result in Fig.<ref> is limited to linear order in $\ell$ or $k$. However, it is important to note that the results depicted in Fig.<ref> and Fig.<ref> hold true for all orders in the expansion of small $\ell$ 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  <cit.> and section-(7.4) of <cit.>. For an un-amputated 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.<ref> and the generalized identity in Fig.<ref>. The identity in Fig.<ref> is a straight forward generalization of the example discussed in Fig.<ref>, which has been tested with the covariantized vertices up to linear order in photon momenta expansion for $\Gamma^{(n)}$ and $\Gamma^{(n+1)}$ vertices with $n=3$.

IR-finite amplitudes:
As we defined earlier, $\mathcal{A}^{(N)}$ represents the all loop scattering amplitude with $N$ number of external massive spinning particles, and $\mathcal{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 $\lbrace e_i\rbrace$ then the following identity holds
𝒬_α_i ^β_i ϵ_iβ_i(p_i)= e_i ϵ_iα_i(p_i) ⇒ 𝒬 ϵ_i(p_i)=e_iϵ_i(p_i) ,  ϵ^T_i(p_i)𝒬^T=e_iϵ^T_i(p_i) .
The K-photon insertion Ward identities of Fig.<ref> and Fig.<ref> imply the exponentiation of the one-loop K-photon contribution $K_{em}$ in (<ref>), as proven in <cit.>. However, the proof of this exponentiation in <cit.> is valid only when the tree-level amplitude of the massive EFT $\left(\Gamma^{(N)}\right)$ is independent of the momenta of scattering particles i.e. it is described by a momenta-independent contact interaction between $N$ number of massive fields, as considered in <cit.>. 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.[We are thankful to P.V. Athira for extensive discussion on this topic.] 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 (<ref>) is
𝒜^(N)≡exp{K_em} 𝒜^(N)_IR-finite  ,

𝒜^(N+1)≡exp{K_em} 𝒜^(N+1)_IR-finite ,
K_em   =      i2∑_i=1^N  ∑_j=1
j≠i^Ne_ie_j ∫d^4ℓ(2π)^4  1ℓ^2-iϵ (2p_i-ℓ)·(2p_j+ℓ) (2p_i·ℓ-ℓ^2+iϵ)
(2p_j·ℓ+ℓ^2-iϵ) .
In equations (<ref>) and (<ref>), $\mathcal{A}^{(N)}_{\text{IR-finite}}$ and $\mathcal{A}^{(N+1)}_{\text{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 (<ref>) and (<ref>). Both $\mathcal{A}^{(N)}_{\text{IR-finite}}$ and $\mathcal{A}^{(N+1)}_{\text{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 K-photon 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
𝒜^(N)_IR-finite,0= Γ^(N) , 𝒜^(N+1)_IR-finite,0=Γ^(N+1) ,
where in the subscript `$0$' corresponds to $0$-loop i.e. tree level. At one-loop order, $\mathcal{A}^{(N)}_{\text{IR-finite}}$ and $\mathcal{A}^{(N+1)}_{\text{IR-finite}}$ are given by:
𝒜^(N)_IR-finite,1   ≡   [𝒜^(N)_G,1+𝒜^(N)_K-finite,1+𝒜^(N)_self,1+𝒜^(N)_non-div,1] ,

𝒜^(N+1)_IR-finite,1   ≡   [𝒜^(N+1)_G,1+𝒜^(N+1)_K-finite,1+𝒜^(N+1)_self,1+𝒜^(N+1)_non-div,1] ,
where in the subscript `$1$' corresponds to $1$-loop. Above

* $\mathcal{A}^{(N)}_{G,1}$ corresponds to the diagram in Fig.<ref> which we need to evaluate with G-photon propagator.

* $\mathcal{A}^{(N)}_{K-\text{finite},1}$ corresponds to the contribution from the diagram in Fig.<ref>, evaluated with K-photon propagator and then subtracted the contribution of $K_{em}\times \Gamma^{(N)}$ at the integrand level.

* $\mathcal{A}^{(N)}_{\text{self,1}}$ represents the set of Feynman diagrams in Fig.<ref> evaluated with full photon propagator.

* $\mathcal{A}^{(N)}_{\text{non-div,1}}$ represents the set of diagrams in Fig.<ref>, evaluated with full photon propagator.

* $\mathcal{A}^{(N+1)}_{G,1}$ corresponds to the set of diagrams in Fig.<ref> when we evaluate them with G-photon propagator.

* $\mathcal{A}^{(N+1)}_{K-\text{finite},1}$ corresponds to the contribution from the diagrams in Fig.<ref>, evaluated with K-photon propagator and then subtracted the contribution $K_{em}\times \Gamma^{(N+1)}$ at the integrand level.

* $\mathcal{A}^{(N+1)}_{\text{self,1}}$ represents the set of Feynman diagrams in Fig.<ref>, evaluated with full photon propagator.

* $\mathcal{A}^{(N+1)}_{\text{non-div,1}}$ represents the set of diagrams in Fig.<ref>, evaluated with full photon propagator.

Diagram contributes to $\mathcal{A}^{(N)}_{G,1}$, where the virtual photon propagator is the G-photon propagator connected between two external hard particle lines. The diagram also contributes to $\mathcal{A}^{(N)}_{K-\text{finite},1}$ when evaluated with K-photon propagator and subtracted $K_{em}\Gamma^{(N)}$ from it's contribution.

§.§ Derivation of soft photon theorem

The goal here will be to derive the order $\omega^{-1}$ and $\ln\omega$ soft factors from the ratio of $\mathcal{A}^{(N+1)}_{\text{IR-finite}}$ and $\mathcal{A}^{(N)}_{\text{IR-finite}}$ when the external photon energy is small i.e. $\omega<<|p_i^\mu |$.

§.§.§ IR-finite one loop $N$-particle amplitude

First we want to analyze all the Feynman diagrams contributing to (<ref>) and show that the result is IR-divergence free. We also provide an explicit integral expression of IR-finite 1-loop amplitude for  $\mathcal{A}^{(N)}_{\text{IR-finite,1}}$.  Using the Feynman rules derived in section-<ref>, the diagram in Fig.<ref> with G-photon propagator contributes to the following
𝒜^(N)_G,1   =   ∑_i,j=1
i> j^N ∫d^4ℓ(2π)^4 -iℓ^2-iϵ 1(p_i-ℓ)^2+m_i^2-iϵ 1(p_j+ℓ)^2+m_j^2 -iϵ

[ ϵ_i^T(-p_i)Γ^(3)_μ(p_i,-p_i+ℓ,-ℓ)Ξ_i(-p_i+ℓ)]  G_(ij)^μν(ℓ)

×[ϵ_j^T(-p_j)Γ^(3)_ν(p_j,-p_j-ℓ,ℓ)Ξ_j(-p_j-ℓ)] Γ^(N)_(ij)(p_i-ℓ,p_j+ℓ)

=    ∑_i,j=1
i> j^N8e_i e_jϵ_i^Tϵ_j^T ∫d^4ℓ(2π)^4 -iℓ^2-iϵ 1(p_i-ℓ)^2+m_i^2-iϵ 1(p_j+ℓ)^2+m_j^2 -iϵ

×{[ p_i^μℓ^ν𝒩^j_μν(-p_j) - p_j^μℓ^ν𝒩^i_μν(-p_i)]  Γ^(N)_(ij)(p_i,p_j)+𝒪(ℓℓ)} ,
where to get the last two lines we have used the identity in (<ref>) for both $i$-th and $j$-th particles and Taylor expanded the numerator in the limit $|\ell^\mu|<<|p_i^\mu| ,|p_j^\mu|$. The expression for $\mathcal{N}^i(-p_i)$ is given in (<ref>). On the other hand the contribution of $\mathcal{A}^{(N)}_{K-\text{finite},1}$ from the diagram in Fig.<ref> becomes
𝒜^(N)_K-finite,1   =    ∑_i,j=1
i> j^N ∫d^4ℓ(2π)^4 -iℓ^2-iϵ 1(p_i-ℓ)^2+m_i^2-iϵ 1(p_j+ℓ)^2+m_j^2 -iϵ

[ ϵ_i^T(-p_i)Γ^(3)_μ(p_i,-p_i+ℓ,-ℓ)Ξ_i(-p_i+ℓ)]  K_(ij)^μν(ℓ)

×[ϵ_j^T(-p_j)Γ^(3)_ν(p_j,-p_j-ℓ,ℓ)Ξ_j(-p_j-ℓ)] Γ^(N)_(ij)(p_i-ℓ,p_j+ℓ)

- K_em×Γ^(N) .
Now evaluating this expression using the identity in (<ref>) and Taylor expanding $\Gamma^{(N)}_{(ij)}(p_i-\ell,p_j+\ell)$ in the limit $|\ell^\mu|<<|p_i^\mu| ,|p_j^\mu|$ we get
𝒜^(N)_K-finite,1   =    i
i> j^Ne_ie_j ϵ_i^Tϵ_j^T∫d^4ℓ(2π)^4  1ℓ^2-iϵ (2p_i-ℓ)·(2p_j+ℓ) (2p_i·ℓ-ℓ^2+iϵ)

×[ -ℓ^ρp_i^ρΓ^(N)_(ij)(p_i,p_j) +ℓ^ρp_j^ρΓ^(N)_(ij)(p_i,p_j) +𝒪(ℓℓ)] .

Diagram contributes to $\mathcal{A}^{(N)}_{\text{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.

Diagrams in Fig.<ref> renormalizes the massive spinning particle propagators in presence of electromagnetic interaction, and all the loops are IR-finite. Say the three diagrams in Fig.<ref> contributes to
ϵ_i^TF_1 Ξ_i(-p_i)p_i^2+m_i^2Γ^(N)_(i)(p_i)  ,  ϵ_i^TF_2 Ξ_i(-p_i)p_i^2+m_i^2Γ^(N)_(i)(p_i)  ,  ϵ_i^TCΞ_i(-p_i)p_i^2+m_i^2Γ^(N)_(i)(p_i)
respectively. Where $F_1,F_2, C$ are unknown constant matrices, which are related using on-shell renormalization condition
C=-F_1 -F_2 .
Hence the on-shell renormalization condition implies
𝒜^(N)_self,1=0 .

Diagram contributes to $\mathcal{A}^{(N)}_{\text{non-div,1}}$, where the virtual photon propagator is the full photon propagator whose one or both end connected to $\Gamma^{(N)}$. In the last diagram the cross corresponds to counter term, which cancels the UV divergences.

In Fig.<ref> 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 $\Gamma^{(N)}$. These diagrams are also IR-finite. For example the first diagram in Fig.<ref> with full photon propagator becomes
∫d^4ℓ(2π)^4 -iℓ^2 -iϵ 1(p_i-ℓ)^2+m_i^2 -iϵ η^μν

×ϵ_i^T(-p_i)[ Γ^(3)_μ(p_i,-p_i+ℓ,-ℓ)Ξ_i(-p_i+ℓ)Γ^(N+1)_(i)ν(p_i-ℓ;ℓ)] .
Now using the identity in (<ref>) for $i$-th particle and the Feynman rule of (<ref>) it is evident that in the limit when the loop momentum $\ell^\mu\rightarrow 0$ the integration behaves as $\int \f{d^4\ell}{|\ell|^3}$, hence is IR-finite. This is the reason we call the set of diagrams in Fig.<ref> as $\mathcal{A}^{(N)}_{\text{non-div,1}}$ as those are IR-divergence free. We do not need to evaluate the contribution $\mathcal{A}^{(N)}_{\text{non-div,1}}$ explicitly for deriving soft photon theorem. Now summing over the contributions of (<ref>), (<ref>), (<ref>) and $\mathcal{A}^{(N)}_{\text{non-div,1}}$, we get
𝒜^(N)_IR-finite,1   =    ∑_i,j=1
i> j^N8e_i e_jϵ_i^Tϵ_j^T ∫d^4ℓ(2π)^4 -iℓ^2-iϵ 1(p_i-ℓ)^2+m_i^2-iϵ 1(p_j+ℓ)^2+m_j^2 -iϵ

×{[ p_i^μℓ^ν𝒩^j_μν(-p_j) - p_j^μℓ^ν𝒩^i_μν(-p_i)]  Γ^(N)_(ij)(p_i,p_j)+𝒪(ℓℓ)}

+i ∑_i,j=1
i> j^Ne_ie_j ϵ_i^Tϵ_j^T∫d^4ℓ(2π)^4  1ℓ^2-iϵ (2p_i-ℓ)·(2p_j+ℓ) (2p_i·ℓ-ℓ^2+iϵ)

×[ -ℓ^ρp_i^ρΓ^(N)_(ij)(p_i,p_j) +ℓ^ρp_j^ρΓ^(N)_(ij)(p_i,p_j) +𝒪(ℓℓ)]

+ 𝒜^(N)_non-div,1 .
From the above expression, it is clear that the loop integrals become infrared finite in the limit as $\ell$ 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).

Set of 1-loop diagrams contributing to $\mathcal{A}^{(N+1)}_{G,1}$, 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}\Gamma^{(N+1)}$ from it's contribution, it also contributes to $\mathcal{A}^{(N+1)}_{K-\text{finite},1}$ .

§.§.§ IR-finite one loop $(N+1)$-particle amplitude in the soft limit

Here we analyze all the Feynman diagrams contributing to (<ref>) in the soft limit i.e. $\omega\rightarrow 0$.
Let us start analyzing the first diagram in Fig.<ref> with G-photon propagator, which has the following expression after using Feynman rules
A_I   ≡    ∑_i=1^N∑_j=1
j≠i^N∫d^4ℓ(2π)^4 -iℓ^2-iϵ 1(p_i+k)^2+m_i^2-iϵ1(p_i+k-ℓ)^2+m_i^2-iϵ1(p_j+ℓ)^2+m_j^2 -iϵ

[ ϵ_i^T(-p_i)ε^ρ(k)Γ^(3)_ρ(p_i,-p_i-k,k)Ξ_i(-p_i-k)Γ^(3)_μ(p_i+k,-p_i-k+ℓ,-ℓ)Ξ_i(-p_i-k+ℓ)]

× G_(ij)^μν(ℓ) [ϵ_j^T(-p_j)Γ^(3)_ν(p_j,-p_j-ℓ,ℓ)Ξ_j(-p_j-ℓ)] Γ^(N)_(ij)(p_i+k-ℓ,p_j+ℓ) .
Substituting the results from (<ref>) and (<ref>) in the above expression and after some manipulation we get
=    ∑_i=1^N∑_j=1
j≠i^Ne_i^2 e_j ϵ_i^Tϵ_j^T1(p_i+k)^2+m_i^2-iϵ∫d^4ℓ(2π)^4 -iℓ^2-iϵ 1(p_i+k-ℓ)^2+m_i^2-iϵ1(p_j+ℓ)^2+m_j^2 -iϵ

[ 8ε.p_i p_j.k +16ε.p_i p_i^ρℓ^σ𝒩^j_ρσ(-p_j)-16ε.p_i p_j^ρℓ^σ𝒩^i_ρσ(-p_i)+8p_i.p_j (ε^ρk^σ-ε^σk^ρ)𝒩^i_ρσ(-p_i)

+4ip_i.k ε^ρp_j^σ𝒦_i(-p_i)p_i^ρΞ_i(-p_i)p_i^σ-4p_i.p_j2p_i.ℓ-ℓ^2+iϵ{4ε.p_i ℓ.k +4p_i.ℓ(ε^ρk^σ-ε^σk^ρ)𝒩^i_ρσ(-p_i)

+2ip_i.kε^ρℓ^σ𝒦_i(-p_i)p_i^ρΞ_i(-p_i)p_i^σ} +𝒪(ℓℓ, kℓ,kk)] {Γ^(N)_(ij)(p_i,p_j)+𝒪(ℓ,k)} .
Inside the square bracket of the numerator in the above expression, we only keep the terms up to linear order in $\ell$ or $k$, as our vertices are derived only up to that order. Note that the above expression is IR-finite in the limit $\ell^\mu\rightarrow 0$ and in the region of integration $|\ell^\mu|<<\omega<<|p_i^\mu |,|p_j^\mu|$ it contributes at order $\mathcal{O}(\omega^0)$. Now to extract $\ln\omega$ contribution we approximate the integrand in the integration range $\omega<<|\ell^\mu|<<|p_i^\mu |,|p_j^\mu|$, by doing so we can approximate
1 (p_i+k-ℓ)^2+m_i^2-iϵ≃1 (p_i-ℓ)^2+m_i^2-iϵ[1+p_i.kp_i.ℓ+iϵ] .
With this approximation the order $\omega^{-1}$ and $\ln\omega$ contribution turns out to be,
=    ∑_i=1^N∑_j=1
j≠i^Ne_iε.p_ip_i.k 8e_i e_j ϵ_i^Tϵ_j^T∫_reg  d^4ℓ(2π)^4 -iℓ^2-iϵ 1(p_i-ℓ)^2+m_i^2-iϵ

1(p_j+ℓ)^2+m_j^2 -iϵ [p_i^ρℓ^σ𝒩^j_ρσ(-p_j)- p_j^ρℓ^σ𝒩^i_ρσ(-p_i)]Γ^(N)_(ij)(p_i,p_j)

j≠i^Ne_i^2 e_j 1p_i.kϵ_i^Tϵ_j^T∫_reg d^4ℓ(2π)^4 1ℓ^2-iϵ 1p_i.ℓ+iϵ1p_j.ℓ-iϵ

[ ε.p_i p_j.k +p_i.kp_i.ℓ+iϵ{2ε.p_i p_i^ρℓ^σ𝒩^j_ρσ(-p_j)-2ε.p_i p_j^ρℓ^σ𝒩^i_ρσ(-p_i)}

+i2p_i.k ε^ρp_j^σ𝒦_i(-p_i)p_i^ρΞ_i(-p_i)p_i^σ-p_i.p_jp_i.ℓ+iϵ{ε.p_i ℓ.k

+i2p_i.kε^ρℓ^σ𝒦_i(-p_i)p_i^ρΞ_i(-p_i)p_i^σ} ] Γ^(N)_(ij)(p_i,p_j) + 𝒪(ω^0).
Above the subscript “$reg$” in the loop integral corresponds to the restricted loop-momentum range $\omega<<|\ell^\mu|<<|p_i^\mu |,|p_j^\mu|$.

Using Feynman rules the second diagram in Fig.<ref> with G-photon propagator becomes
A_II   ≡    ∑_i=1^N∑_j=1
j≠i^N∫d^4ℓ(2π)^4 -iℓ^2-iϵ 1(p_i+k-ℓ)^2+m_i^2-iϵ1(p_j+ℓ)^2+m_j^2 -iϵ

[ ϵ_i^T(-p_i)ε^ρ(k)Γ^(4)_ρμ (p_i,-p_i-k+ℓ,k,-ℓ)Ξ_i(-p_i-k+ℓ)] G_(ij)^μν(ℓ)

×  [ϵ_j^T(-p_j)Γ^(3)_ν(p_j,-p_j-ℓ,ℓ)Ξ_j(-p_j-ℓ)] Γ^(N)_(ij)(p_i+k-ℓ,p_j+ℓ).
After substituting the results from (<ref>) and (<ref>) and keeping terms which can contribute up to order $\ln\omega$ in the integration range $reg\equiv (\omega<<|\ell^\mu|<<|p_i^\mu|, |p_j^\mu|)$ we get,
A_II   =   -i∑_i=1^N∑_j=1
j≠i^Ne_i^2e_j ϵ_i^Tϵ_j^T∫_reg d^4ℓ(2π)^4 1ℓ^2-iϵ 1p_i.ℓ+iϵ1p_j.ℓ-iϵ

( ε.ℓ+i2ε^ρℓ^σ{𝒦_i(-p_i)p_i^ρΞ_i(-p_i)p_i^σ+𝒦_i(-p_i)p_i^σΞ_i(-p_i)p_i^ρ})] Γ^(N)_(ij)(p_i,p_j)

+ 𝒪(ω^0) .

Set of 1-loop diagrams contributing to $\mathcal{A}^{(N+1)}_{\text{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.<ref> with G-photon propagator becomes
A_III   ≡    ∑_i=1^N∑_j=1
j≠i^N∫d^4ℓ(2π)^4 -iℓ^2-iϵ 1(p_i-ℓ)^2+m_i^2-iϵ1(p_i+k-ℓ)^2+m_i^2-iϵ1(p_j+ℓ)^2+m_j^2 -iϵ

[ ϵ_i^T(-p_i)ε^ρ(k)Γ^(3)_μ(p_i,-p_i+ℓ,-ℓ)Ξ_i(-p_i+ℓ)Γ^(3)_ρ(p_i-ℓ,-p_i-k+ℓ,k)Ξ_i(-p_i-k+ℓ)]

× G_(ij)^μν(ℓ) [ϵ_j^T(-p_j)Γ^(3)_ν(p_j,-p_j-ℓ,ℓ)Ξ_j(-p_j-ℓ)] Γ^(N)_(ij)(p_i+k-ℓ,p_j+ℓ) .
After substituting the results from (<ref>) and (<ref>) and keeping terms which can contribute up to order $\ln\omega$ in the integration range $\omega<<|\ell^\mu|<<|p_i^\mu|,|p_j^\mu|$ we get,
A_III   =    -i∑_i=1^N∑_j=1
j≠i^Ne_i^2e_j ϵ_i^Tϵ_j^T∫_reg d^4ℓ(2π)^4 1ℓ^2-iϵ 1(p_i.ℓ+iϵ)^21p_j.ℓ-iϵ

[ 2ε.p_i p_i^ρℓ^σ𝒩^j_ρσ(-p_j)-2ε.p_i p_j^ρℓ^σ𝒩^i_ρσ(-p_i)-i2p_i.ℓp_j^ρε^σ𝒦_i(-p_i)p_i^ρΞ_i(-p_i)p_i^σ

+i2p_i.p_j ℓ^ρε^σ𝒦_i(-p_i)p_i^ρΞ_i(-p_i)p_i^σ] Γ^(N)_(ij)(p_i,p_j) + 𝒪(ω^0) .
The fourth diagram in Fig.<ref> with G-graviton contributes at order $\omega^{-1}$ but won't contribute at order $\ln\omega$. The order $\omega^{-1}$ contribution turns out to be,
=    ∑_k=1^Ne_kε.p_kp_k.k ∑_i=1
j>i^N8e_i e_j ϵ_i^Tϵ_j^T∫_reg  d^4ℓ(2π)^4 -iℓ^2-iϵ 1(p_i-ℓ)^2+m_i^2-iϵ

1(p_j+ℓ)^2+m_j^2 -iϵ [p_i^ρℓ^σ𝒩^j_ρσ(-p_j)- p_j^ρℓ^σ𝒩^i_ρσ(-p_i)]Γ^(N)_(ij)(p_i,p_j) + 𝒪(ω^0) .
The fifth diagram in Fig.<ref> start contributing from order $\omega^0$ in the soft expansion when we evaluate it with G-photon propagator i.e. $A_{V}=0+\mathcal{O}(\omega^{0})$. Here we are not writing down the non-vanishing contribution of $A_{V}$ at order $\omega^{0}$ explicitly, as it is not essential for deriving the order $\omega^{-1}$ and $\ln\omega$ soft factors. Now summing over the external particle legs, total contribution of $\mathcal{A}^{(N+1)}_{G,1}$ at orders $\omega^{-1}$ and $\ln\omega$ turns out to be

=    A_I+A_II+A_III+A_IV+A_V

=    ∑_k=1^Ne_kε.p_kp_k.k∑_i=1^N  ∑_j=1
j>  i^N 8e_i e_j ϵ_i^Tϵ_j^T∫_reg d^4ℓ(2π)^4 -iℓ^2-iϵ 1(p_i-ℓ)^2+m_i^2-iϵ

1(p_j+ℓ)^2+m_j^2 -iϵ [p_i^ρℓ^σ𝒩^j_ρσ(-p_j)- p_j^ρℓ^σ𝒩^i_ρσ(-p_i)]Γ^(N)_(ij)(p_i,p_j)

+ i ∑_i=1^N  ∑_j=1
j≠i^Ne_i^2e_j ϵ_i^Tϵ_j^T∫_reg d^4ℓ(2π)^4 1ℓ^2-iϵ 1p_i.ℓ+iϵ 1p_j.ℓ-iϵ

×[ ε.p_ip_i.kp_j.k -ε.p_ip_i.k p_i.p_j ℓ.kp_i.ℓ+iϵ-ε.p_j +ε.ℓp_i.p_jp_i.ℓ+iϵ]Γ^(N)_(ij)(p_i,p_j) +𝒪(ω^0) .
Interestingly in the above expression all the theory dependent pieces involving $\mathcal{K}_i\ , \Xi_i\ , \mathcal{B}^i$ at order $\ln\omega$ cancels out at the integrand level, when we some over the contributions (<ref>),(<ref>),(<ref>). This confirms the fact that the subleading soft photon theorem at order $\ln\omega$ is universal. The above result can also be rewritten in the following compact form
𝒜^(N+1)_G,1    =   ∑_k=1^Ne_kε.p_kp_k.k  𝒜^(N)_G,1

+   ∑_i=1^Ne_iε_μk_νp_i·k{(p_i^μp_iν-p_i^νp_iμ)K_em^reg} Γ^(N) +𝒪(ω^0) ,
where the expression of $\mathcal{A}^{(N)}_{G,1}$ is given in (<ref>), and $K_{em}^{reg}$ is the approximated form of the integral $K_{em}$ in (<ref>) in the integration range $\omega<<|\ell^\mu|<<|p_i^\mu|,|p_j^\mu|$. The integration has been explicitly evaluated in <cit.> and the result reads
K_em^reg≡i2∑_ℓ=1^N  ∑_j=1
j≠ℓ^Ne_ℓe_j ∫_reg d^4ℓ(2π)^4  1ℓ^2-iϵ (p_ℓ·p_j)(p_ℓ·ℓ+iϵ) (p_j·ℓ- iϵ)

≃    -i2∑_ℓ=1^N  ∑_j=1
j≠ℓ^Ne_ℓe_j4π (lnω)  p_ℓ·p_j√((p_ℓ.p_j)^2-p_ℓ^2p_j^2) {δ_η_ℓη_j,1-i2πln(p_ℓ.p_j+√((p_ℓ.p_j)^2-p_ℓ^2p_j^2)p_ℓ.p_j-√((p_ℓ.p_j)^2-p_ℓ^2p_j^2))} .

Above $\eta_j=+1$ if $j$-th particle is ingoing and $\eta_j=-1$ if $j$-th particle is outgoing and under $\simeq$ sign we only have written the order $\ln\omega$ contribution while evaluating the integral.

The first diagram above corresponds to the sum over the contributions of the first four diagrams in Fig.<ref> when we replace the photon polarization $\varepsilon^\mu$ by $k^\mu$. The second diagram above corresponds to the sum over the contributions of the fifth and sixth diagrams in Fig.<ref> when we replace the photon polarization $\varepsilon^\mu$ by $k^\mu$.

To evaluate the contribution of $\mathcal{A}^{(N+1)}_{K-\text{finite},1}$, first we compute the diagrams in Fig.<ref> with the virtual photon propagator replaced by the K-photon propagator. Afterward, we subtract the contribution of $K_{em}\Gamma^{(N+1)}$ from this computed result. The final contribution of $\mathcal{A}^{(N+1)}_{K-\text{finite},1}$ within the regulated range of integration $\omega<<|\ell^\mu|<<|p_i^\mu|,|p_j^\mu|$ becomes
𝒜^(N+1)_K-finite,1   =     ∑_k=1^Ne_kε.p_kp_k.k 𝒜^(N)_K-finite,1 +𝒪(ω^0) ,
where the expression of $\mathcal{A}^{(N)}_{K-\text{finite},1}$ is given in (<ref>). Note that the above result contributes at order $\omega^{-1}$ and does not contribute to order $\ln\omega$ soft theorem.

The self-energy kind of diagrams contributing to $\mathcal{A}^{(N+1)}_{\text{self,1}}$ in Fig.<ref> are not necessary to compute explicitly as these diagrams sum up to zero using on-shell renormalization condition (<ref>) as we are going to discuss below. The contribution from the sum of the first four diagrams in Fig.<ref> can be described by the following general structure $\epsilon_i^{T}\varepsilon.p_i f_{1}(p_i.k)\Gamma_{(i)}^{(N)}(p_i+k)$. Similarly the contribution from the sum of fifth and sixth diagrams in Fig.<ref> can be described by the following structure $\epsilon_i^{T}\varepsilon.p_i f_{2}(p_i.k)\Gamma_{(i)}^{(N)}(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 $\varepsilon^\mu \rightarrow k^\mu$ and using Ward identity. Using the diagrammatic identities of Fig.<ref> and Fig.<ref>,  the first four diagrams in Fig.<ref> after replacing $\varepsilon^\mu \rightarrow k^\mu$ reduces to the first diagram in Fig.<ref>. Similarly using the diagrammatic identities of Fig.<ref> and Fig.<ref>,  the sum of fifth and sixth diagrams in Fig.<ref> after replacing $\varepsilon^\mu \rightarrow k^\mu$ reduces to the second diagram in Fig.<ref>. Now using the constant matrices $F_1,F_2$ introduced in (<ref>) for the diagrams in Fig.<ref> and comparing the general structure above we get
f_1(p_i.k)=-12(p_i.k)^2 F_1𝒬_iΞ_i(-p_i-k) , f_1(p_i.k)=-12(p_i.k)^2 F_2𝒬_iΞ_i(-p_i-k) .
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.<ref> we get
ϵ_i^Tε.p_i2(p_i.k)^2[-F_1𝒬_i -F_2𝒬_i +𝒬_i^TC] Ξ_i(-p_i-k)Γ^(N)_(i)(p_i+k) .
Note that $F_1$, $F_2$ satisfy the same property under the operation of charge matrix $\mathcal{Q}$  as $\Xi_i$ satisfies in (<ref>). Hence using this relation the terms inside the square bracket in the above expression reduces to $\mathcal{Q}_i^{T}(F_1+F_2+C)$ which vanishes using the on-shell renormalization condition (<ref>). Note that in the above expression we neglected the possible order $\omega^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.<ref> also vanishes up to possible order $\omega^0$ contribution . Hence in the soft limit, the sum over all the Feynman diagrams in Fig.<ref> contribute as
𝒜^(N+1)_self,1  = 0 +𝒪(ω^0) .

Set of 1-loop diagrams contributing to $\mathcal{A}^{(N+1)}_{\text{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.

The diagrams in Fig.<ref> contributing to $\mathcal{A}^{(N+1)}_{\text{non-div,1}}$ are IR-divergence free, as in $\ell^\mu\rightarrow 0$ limit and finite $k^\mu$ the third, sixth and seventh diagrams behaves like $\int \f{d^4\ell}{|\ell|^3}$, and first, second, fourth and fifth diagrams behaves like $\int \f{d^4\ell}{|\ell|^2}$. On the other hand the sum of first, fourth and sixth diagrams in Fig.<ref> contribute to leading soft factor at order $\omega^{-1}$. Now in the integration region $\omega<<|\ell^\mu|<<|p_i^\mu |$, individually the first and third diagrams in Fig.<ref> behave like $\int_{reg} \f{d^4\ell}{\ell^2-i\epsilon}\f{1}{(p_i.\ell+i\epsilon )^2} $ after the expansion of the propagators, hence those have the possibility of contributing at order $\ln\omega$. But when we sum over the contributions of the first, second and third diagrams, the order $\ln\omega$ contributing coefficient of $\int_{reg} \f{d^4\ell}{\ell^2-i\epsilon}\f{1}{(p_i.\ell+i\epsilon )^2} $ cancels each other, and  left out part starts contributing from order $\omega^0$. Hence leaving first, fourth and sixth diagrams, all the other diagrams start contributing at order $\omega^0$ in the soft expansion. After summing over all the contributions, we get
𝒜^(N+1)_non-div,1=∑_i=1^N e_iε.p_ip_i.k 𝒜^(N)_non-div,1 +𝒪(ω^0) .

Soft photon theorem result at one-loop:
Summing over the contributions of

(<ref>),(<ref>),(<ref>),(<ref>) in the soft limit, we get the following soft theorem expression
𝒜^(N+1)_IR-finite,1   =   [𝒜^(N+1)_G,1+ 𝒜^(N+1)_K-finite,1+𝒜^(N+1)_self,1+𝒜^(N+1)_non-div,1]

ω→0 =     ∑_i=1^Ne_iε.p_ip_i.k  𝒜^(N)_IR-finite,1

+   ∑_i=1^Ne_iε_μk_νp_i·k{(p_i^μp_iν-p_i^νp_iμ)K_em^reg} 𝒜^(N)_IR-finite,0 +𝒪(ω^0) .
This result agrees with the loop corrected subleading soft photon theorem, originally derived in <cit.> for minimally coupled scalar QED. The derivation of this result further confirms the universality of the $\ln\omega$ 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 $\omega\ln\omega$ soft factor from the aforementioned one-loop amplitude, we require the vertices: $\Gamma^{(3)}$ as given in equation (<ref>) up to order $\mathcal{O}(\ell^2)$, $\Gamma^{(4)}$ as given in equation (<ref>) up to order $\mathcal{O}(\ell_1, \ell_2)$, and $\widetilde{\Gamma}^{(N+1)}$ up to order $\mathcal{O}(\ell)$. However, when dealing with non-minimal couplings as expressed in equation (<ref>), it is currently unknown how to derive $\widetilde{\Gamma}^{(N+1)}$ in terms of $\Gamma^{(N)}$ at order $\mathcal{O}(\ell)$. Additionally, at order $\mathcal{O}(\ell^2)$, new sets of non-minimal couplings will contribute to $\Gamma^{(3)}$. Due to these reasons, we are unable to derive the order $\omega\ln\omega$ soft photon theorem in this article, and it is not clear whether such a soft factorization at order $\omega\ln\omega$ is possible or not.

In the theory of scalar-QED minimally coupled to gravity the order $\ln\omega$ correction to (<ref>) due to gravitational interaction has also been derived in <cit.> by analyzing one-loop amplitudes. In the theory of scalar-QED soft photon theorem at order $\omega (\ln\omega)^2$ has also been derived in <cit.> analyzing two-loop amplitudes and the soft factor is provided in (<ref>). Multiple soft photon theorem up to subleading order in soft expansion has also been derived in section-(3.5) of <cit.>, by analyzing one-loop amplitudes and the result reads:

ω→0 =  ∏_ℓ=1^M{∑_i=1^Ne_iε_ℓ·p_ip_i·k_ℓ} 𝒜^(N)_IR-finite,1

m≠ℓ^M{∑_j=1^Ne_jε_m·p_jp_j·k_m}] ∑_i=1^Ne_iε_ℓμ k_ℓνp_i·k_ℓ{(p_i^μp_iν-p_i^νp_iμ)K_em^reg} 𝒜^(N)_IR-finite,0

+ 𝒪(ω^-M+1) .
Above $\mathcal{A}^{(N+M)}_{\text{IR-finite,1}}$ represents the IR-finite part of 1-loop scattering amplitude involving $N$-number of hard particles and $M$-number of soft photons with polarizations and momenta $\lbrace \varepsilon_\ell ,k_\ell=-\omega \mathbf{n}_\ell\rbrace$ for $\ell=1,\cdots, M$ as external states.

§ 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 <cit.> for a generic theory of quantum gravity. This derivation will establish the universal (theory independent) nature of  $\ln\omega$ soft graviton factor.

§.§ Sen's Covariantization prescription and Feynman rules

The covariantization of the quadratic part of the massive EFT action (<ref>) in the soft gravitational background has been carried out in <cit.>. 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 <cit.>, as described below. In this work, we use a different parametrization of gravitational fluctuation, which is closely related to the one presented in <cit.>. 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
h_μν(x)  ≡  12κ(g_μν(x)-η_μν) ,
where $\kappa=\sqrt{8\pi G}\ $ with $G$ being the four dimensional Newton's constant. In de Donder gauge $\p^\mu h_{\mu\nu}=\f{1}{2}\p_\nu h^{\rho}_\rho$ the graviton propagator becomes
Δ_F^μν,ρσ(ℓ)=-iℓ^2 -iϵ 12(η^μρη^νσ+η^μση^νρ-η^μνη^ρσ) ,
where $(\mu\nu)$ and $(\rho\sigma)$ 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 $\varepsilon_{\mu\nu}(k)$ which satisfies the traceless and transverse conditions
η^μνε_μν(k)=0 , k^με_μν(k)=0 .
With the definition of gravitational fluctuation (<ref>), we express the following quantities as a power series expansion of $\kappa$:
√(-det g)   =   1+κh+κ^22(h^2 -2h^μνh_μν)+κ^36(h^3+8h^μρ h_ρν h^ν_ μ-6hh^μνh_μν)+𝒪(κ^4) ,

g^μν   =   η^μν-2κh^μν+4κ^2 h^μρh_ρ^ ν -8κ^3 h^μτh_τρh^ρν+𝒪(κ^4) ,

e_μ^ a    =   δ_μ^a +κh_μ^ a-κ^22h_ρ^ ah_μ^ ρ+κ^32h_ρ^ ah_ν^ ρh_μ^ ν+𝒪(κ^4) ,

E_a^ μ    =   δ_a^μ-κh_a^ μ+32κ^2 h_a^ ρh_ρ^ μ-52κ^3 h_a^ ρh_ρνh^νμ+𝒪(κ^4) ,

Γ^λ_μν   =   κ(η^λσ-2κh^λσ+4κ^2 h^λρh_ρ^ σ)( _μh_νσ+_νh_σμ-_σh_μν)+𝒪(κ^4) ,

ω_μ^ab   ≡    η^bc e_ν^ aD_μE_c^ ν=η^bc e_ν^ a_μE_c^ ν+η^bc e_ν^ aΓ_μρ^νE_c^ ρ

=    κ(^b h_μ^ a-^a h_μ^ b)+κ^2( 12h^bν_μh_ν^ a-12h^aν_μh_ν^ b+h^bν^a h_μν-h^aν^b h_μν

+h^aν_νh_μ^ b-h^bν_νh_μ^ a)+ 𝒪(κ^3) ,

R_μνρσ   =    κ(_ρ_νh_μσ+_μ_σh_νρ-_μ_ρh_νσ-_σ_νh_μρ)+𝒪(κ^2) .
In the RHS of all the above expressions both curved space indices ($\mu, \nu, \rho, \sigma, \ldots$) and tangent space indices ($a, b, c, \ldots$) are raised or lowered by using the Minkowski metric $\eta$. The trace of the metric fluctuation is defined as $h\equiv h_{\mu\nu}\eta^{\mu \nu}$. Additionally, the symbols $e_\mu^{\ a}$ represent the vierbein, $E_a^{\ \mu}$ represents the inverse vierbein, $\omega_\mu^{ab}$ represents the spin-connection, $\Gamma^{\lambda}_{\mu\nu}$ represents the Christoffel connection, and $R_{\mu\nu\rho\sigma}$ represents the Riemann tensor.

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 $\Phi_\beta$ should be replaced by covariant derivative in the curved space after multiplication of inverse vierbeins in the following way
_a_1_a_2⋯_a_nΦ_β → E_a_1^ μ_1E_a_2^ μ_2⋯E_a_n^ μ_n D_(μ_1D_μ_2⋯D_μ_n)Φ_β ,
where the expression of one, two and three covariant derivatives on $\Phi_\beta$ are given by
D_μΦ_β   =   _μΦ_β+[12ω_μ^ab(Σ_ab)_β^ γΦ_γ] ,

D_(μ D_ν)Φ_β   =   _μ_νΦ_β+12(ω_μ^ab(Σ_ab)_β^ γ_νΦ_γ+ω_ν^ab(Σ_ab)_β^ γ_μΦ_γ)+[12_(μω_ν)^ab(Σ_ab)_β^ γΦ_γ]

+14ω_(μ^abω_ν)^cd(Σ_ab)_β^ γ(Σ_cd)_γ^ δΦ_δ-[Γ^ρ_μν_ρΦ_β] -12Γ^ρ_μνω_ρ^ab(Σ_ab)_β^ γΦ_γ ,

D_μD_νD_ρΦ_β   =    _μ_ν_ρΦ_β+12_μω_ν^ab(Σ_ab)_β^ γD_ρΦ_γ+12 ω_ν^ab(Σ_ab)_β^ γ_μD_ρΦ_γ

-[ _μΓ^σ_νρ_σΦ_β]-12 _μΓ^σ_νρ ω_σ^ab(Σ_ab)_β^ γΦ_γ- Γ^σ_νρ_μD_σΦ_β

+12ω_μ^ab(Σ_ab)_β^ γD_νD_ρΦ_γ-Γ^σ_μνD_σD_ρΦ_β-Γ^σ_μρD_νD_σΦ_β .
The last equation above needs to symmetrize in $\mu,\nu,\rho$ indices before substituting in (<ref>). 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 $\Phi$. Also while covariantizing the action (<ref>), we need to include $\sqrt{-\text{det }g}$ as a covariant measure of volume.

The Fourier transform of gravitational fluctuation $h_{\mu\nu}(x)$ will be denoted by $h_{\mu\nu}(\ell)$ and the relation between them is given by
h_μν(x)≡∫d^4ℓ(2π)^4  e^iℓ.x h_μν(ℓ) .

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
S^(3)   =    κ2∫d^4 q_1(2π)^4d^4 q_2(2π)^4d^4ℓ(2π)^4 (2π)^4δ^(4)(q_1+q_2+ℓ)Φ_α(q_1)[ h(ℓ)𝒦^αβ(q_2)-h_a^μ(ℓ) q_2μ𝒦^αβ(q_2)q_2a

+12 {ℓ^b h^a_μ(ℓ)-ℓ^a h_μ^b(ℓ)}𝒦^αγ(q_2)q_2μ(Σ_ab)_γ^ β+14 ℓ_μ{ℓ^b h^a_ν(ℓ)-ℓ^a h_ν^b(ℓ)} ^2 𝒦^αγ(q_2)q_2μ q_2ν(Σ_ab)_γ^ β

-12{ℓ_μh^ρ_ν(ℓ) +ℓ_νh^ρ_μ(ℓ) -ℓ^ρh_μν(ℓ)}q_2ρ^2 𝒦^αβ(q_2)q_2μ q_2ν

-16ℓ_μ{ℓ_νh_ρ^σ(ℓ) +ℓ_ρh_ν^σ(ℓ) -ℓ^σh_νρ(ℓ)}q_2σ^3 𝒦^αβ(q_2)q_2μ q_2ν q_2ρ +𝒪(ℓ^3)]Φ_β(q_2) .
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
S^(3)= 12       ∫d^4q_1(2π)^4d^4q_2(2π)^4d^4ℓ(2π)^4 (2π)^4δ^(4)(q_1+q_2+ℓ)

Φ_α(q_1) [R_μνρσ(ℓ)  𝒢^αβ,μνρσ(q_2) +𝒪(ℓ^3)] Φ_β(q_2) ,
where $R_{\mu\nu\rho\sigma}(\ell)$ is the Fourier transform of the linearized Riemann tensor in (<ref>) which reads
R_μνρσ(ℓ)=-κ(ℓ_ρℓ_νh_μσ+ℓ_μℓ_σh_νρ-ℓ_μℓ_ρh_νσ-ℓ_σℓ_νh_μρ) .
In the non-minimal action (<ref>), $\mathcal{G}$ satisfies the following property
𝒢^αβ,μνρσ(q_2)   =   ± 𝒢^βα,μνρσ(-q_1-ℓ) .
In the above equation $+$ sign is for Grassmannian even field and $-$ sign is for Grassmannian odd field. Again during the derivation we consider $\Phi$ 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
S^(4)   =   κ^22∫d^4 q_1(2π)^4d^4 q_2(2π)^4d^4ℓ_1(2π)^4d^4ℓ_2(2π)^4 (2π)^4δ^(4)(q_1+q_2+ℓ_1+ℓ_2)

Φ_α(q_1)[ 12{h(ℓ_1)h(ℓ_2)-2h^μν(ℓ_1)h_μν(ℓ_2)}𝒦^αβ(q_2)+h(ℓ_1)(-h_a^μ(ℓ_2) q_2μ𝒦^αβ(q_2)q_2a

+12 {ℓ_2^b h^a_μ(ℓ_2)-ℓ_2^a h_μ^b(ℓ_2)}𝒦^αγ(q_2)q_2μ(Σ_ab)_γ^ β-12{ℓ_2μ h^ρ_ν(ℓ_2) +ℓ_2ν h^ρ_μ(ℓ_2) -ℓ_2^ρh_μν(ℓ_2)}

×q_2ρ ^2 𝒦^αβ(q_2)q_2μ q_2ν)+32h_a^ρ(ℓ_1)h_ρ^μ(ℓ_2)q_2μ𝒦^αβ(q_2)q_2a +12h_a^μ(ℓ_1)h_b^ν(ℓ_2)q_2μq_2ν^2 𝒦^αβ(q_2)q_2a q_2b

-12h_c^μ(ℓ_1){ℓ_2^b h_μ^a(ℓ_2)-ℓ_2^a h_μ^b(ℓ_2)}𝒦^αγ(q_2)q_2c(Σ_ab)_γ^ β-14h_c^μ(ℓ_1)q_2μ{ℓ_2^b h_ν^a(ℓ_2)-ℓ_2^a h_ν^b(ℓ_2)}

×^2 𝒦^αγ(q_2)q_2c q_2ν(Σ_ab)_γ^ β +12h_a^μ(ℓ_1){ℓ_2μh_ν^ρ(ℓ_2)+ℓ_2νh_μ^ρ(ℓ_2)-ℓ_2^ρh_μν(ℓ_2)}q_2ρ^2 𝒦^αβ(q_2)q_2a q_2ν

+16h_a^μ(ℓ_1)q_2μ{ℓ_2νh_σ^ρ(ℓ_2)+ℓ_2σh_ν^ρ(ℓ_2)-ℓ_2^ρh_νσ(ℓ_2)}q_2ρ ^3 𝒦^αβ(q_2)q_2aq_2ν q_2σ

+12{12h^bσ(ℓ_1)ℓ_2μh_σ^a(ℓ_2)-12h^aσ(ℓ_1)ℓ_2μh_σ^b(ℓ_2)+h^bσ(ℓ_1)ℓ_2^a h_μσ(ℓ_2)-h^aσ(ℓ_1)ℓ_2^b h_μσ(ℓ_2)

+h^aσ(ℓ_1)ℓ_2σh_μ^b(ℓ_2)-h^bσ(ℓ_1)ℓ_2σh_μ^a(ℓ_2)}𝒦^αγ(q_2)q_2μ(Σ_ab)_γ^ β +h^ρσ(ℓ_1){ℓ_2μh_νσ(ℓ_2)

+ℓ_2νh_μσ(ℓ_2)-ℓ_2σh_μν(ℓ_2)}q_2ρ^2 𝒦^αβ(q_2)q_2μq_2ν +𝒪(ℓ_1^2, ℓ_1ℓ_2,ℓ_2^2)]Φ_β(q_2) .

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:
S_EFT   =    12κ^2∫d^4x√(-det(g_μν)) (R +𝒪(R^2)) .

Feynman rules for the vertices:
Starting from the interacting parts of the action given in (<ref>) and (<ref>), Feynman rule for the vertex describing interaction between two spinning particles with momenta $q$ and $-(q+\ell)$, and one graviton with Lorentz indices $\mu\nu$ and momentum $\ell$ turns out to be

=    iκ[ η_μν𝒦(-q)+12η_μνℓ^ρ𝒦(-q)q^ρ+14η_μνℓ^ρℓ^σ^2 𝒦(-q)q^ρq^σ- q_(μ 𝒦(-q)q^ν)-12ℓ_(μ 𝒦(-q)q^ν)

-12q_(μℓ^ρ^2 𝒦(-q)q^ν) q^ρ -14q_(μℓ^ρℓ^σ^3 𝒦(-q)q^ν) q^ρq^σ-12ℓ_(μ ℓ^ρ^2 𝒦(-q)q^ν) q^ρ-12ℓ^b𝒦(-q)q^(μΣ_ν)b

+12ℓ^b  Σ^T_(νb 𝒦(-q)q^μ) -14ℓ^b ℓ^ρ ^2 𝒦(-q)q^(μ q^ρΣ_ν) b +14ℓ^b ℓ^ρΣ^T_(νb^2 𝒦(-q)q^μ) q^ρ

+14{ℓ_μ ℓ^ρ^2 𝒦(-q)q^ρq^ν+ℓ_ν ℓ^ρ^2 𝒦(-q)q^ρq^μ-ℓ^2   ^2 𝒦(-q)q^μq^ν}

+112ℓ^σ{q_μ ℓ^ρ^3 𝒦(-q)q^ρq^ν q^σ+q_ν ℓ^ρ^3 𝒦(-q)q^ρq^μ q^σ-q.ℓ^3 𝒦(-q)q^μq^νq^σ}

-ℓ^ρℓ^σ{𝒢_(μρσν)(-q)+𝒢_σ(νμ)ρ(-q)-𝒢_(μρν)σ(-q) -𝒢_σ(νρμ)(-q)}+𝒪(ℓ^3)] ,
where we suppressed the massive particle spin/polarization indices. The above expression is symmetrized under $\mu \leftrightarrow\nu$ 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 (<ref>), Feynman rule for the vertex describing interaction between two spinning particles with momenta $q$ and $-(q+\ell_1+\ell_2)$, and two gravitons with Lorentz indices $(\mu\nu), (\rho\sigma)$ and momenta $\ell_1, \ell_2$ respectively, turns out to be

=   iκ^2 [ (η_μνη_ρσ-2η_μρη_νσ){𝒦(-q)+12(ℓ_1+ℓ_2)^κ𝒦(-q)q^κ}-η_μν{q_ρ𝒦(-q)q^σ

+12(ℓ_1+ℓ_2)_ρ𝒦(-q)q^σ+12q_ρ(ℓ_1+ℓ_2)^κ^2 𝒦(-q)q^σq^κ+12ℓ_2^b 𝒦(-q)q^ρΣ_σb -12ℓ_2^b Σ^T_σb𝒦(-q)q^ρ}

-η_ρσ{q_μ𝒦(-q)q^ν+12(ℓ_1+ℓ_2)_μ𝒦(-q)q^ν+12q_μ(ℓ_1+ℓ_2)^κ^2 𝒦(-q)q^νq^κ+12ℓ_1^b 𝒦(-q)q^μΣ_νb

-12ℓ_1^b Σ^T_νb𝒦(-q)q^μ}+32η_μρ{q_σ𝒦(-q)q^ν+12(ℓ_1+ℓ_2)_σ𝒦(-q)q^ν+12q_σ(ℓ_1+ℓ_2)^κ^2 𝒦(-q)q^κq^ν

+q_ν𝒦(-q)q^σ +12(ℓ_1+ℓ_2)_ν𝒦(-q)q^σ+12q_ν(ℓ_1+ℓ_2)^κ^2 𝒦(-q)q^κq^σ}+q_μq_ρ^2 𝒦(-q)q^νq^σ

+12{q_μ(ℓ_1 +ℓ_2)_ρ+q_ρ(ℓ_1+ℓ_2)_μ}^2 𝒦(-q)q^νq^σ+12q_μq_ρ(ℓ_1 +ℓ_2)^κ^3 𝒦(-q)q^νq^σq^κ

+12η_μρ{ℓ_2^b 𝒦(-q)q^νΣ_σb-ℓ_2^b Σ^T_σb𝒦(-q)q^ν +ℓ_1^b 𝒦(-q)q^σΣ_νb -ℓ_1^b Σ^T_νb𝒦(-q)q^σ}

+14q_μℓ_2^b {^2 𝒦(-q)q^νq^ρΣ_σb -Σ_σb^T ^2 𝒦(-q)q^νq^ρ}+14q_ρℓ_1^b {^2 𝒦(-q)q^σq^μΣ_νb-Σ_νb^T ^2 𝒦(-q)q^σq^μ}

-14η_νσ(ℓ_2-ℓ_1)^κ{𝒦(-q)q^κΣ_ρμ-Σ_ρμ^T𝒦(-q)q^κ}-12η_νσℓ_2^a {𝒦(-q)q^ρΣ_aμ -Σ_aμ^T 𝒦(-q)q^ρ}

-12η_νσℓ_1^a {𝒦(-q)q^μΣ_aρ -Σ_aρ^T 𝒦(-q)q^μ}-12ℓ_2ν{𝒦(-q)q^σΣ_μρ-Σ_μρ^T 𝒦(-q)q^σ}

-12ℓ_1σ{𝒦(-q)q^νΣ_ρμ-Σ_ρμ^T 𝒦(-q)q^ν}+𝒪(ℓ_1^2, ℓ_1ℓ_2,ℓ_2^2)] ,
where we suppressed the massive particle spin/polarization indices. We could have symmetrize the above expression under $\mu\leftrightarrow\nu$ and $\rho\leftrightarrow\sigma$ 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.

Now we want to evaluate the amputated Green's function involving $N$ number of massive spinning particles and one graviton with momentum $\ell$ and Lorentz indices $\mu\nu$, where the graviton is not connected to any external particle leg. Using covariantization prescription it has been evaluated in appendix-<ref>, where additional complication has been taken care off due to the presence of different momentum conserving delta functions between $\widetilde{\Gamma}^{(N+1)}$ and $\Gamma^{(N)}$. The final result reads from (<ref>) turns out to be
Γ^(N+1)α_1⋯α_N_μν(ℓ)   =    κ∑_i=1^N[-δ_β_i^α_i p_i(μp_i^ν) +(Σ_i(νb)_β_i^ α_iℓ^b  p_i^μ) -12δ_β_i^α_i{ℓ^ρp_iμ^2p_i^ρp_i^ν

The three graviton vertex follows from (<ref>) with the ingoing graviton momenta $k$, $\ell $ and $-(\ell+k)$, and the Lorentz indices carried by them $(\mu\nu)$, $(\rho\sigma)$ and $(\tau\kappa)$ respectively takes the following form:
V^(3)_μν, ρσ,τκ(k,ℓ, -ℓ-k)

=    i  κ [  -(k.ℓ+ℓ^2+k^2)( η_μνη_ρτη_σκ+

+ 4(k.ℓ+ℓ^2+k^2)η_ρνη_στη_κμ+12(k.ℓ+ℓ^2 +k^2)η_μνη_ρση_τκ

- 2  (k_τℓ_κη_μρη_νσ-ℓ_μ(k+ℓ)_νη_ρτη_σκ-(k+ℓ)_ρ k_ση_μτη_νκ)

k_τ +ℓ_ρk_τ-ℓ_τk_ρ)η_μση_νκ)] +𝒪(k^4, k^3ℓ, k^2ℓ^2, kℓ^3,ℓ^4) .

§.§ KG-decomposition and IR-finite amplitudes

In Grammer-Yennie prescription we decompose the internal graviton propagator (<ref>) with momentum $\ell$ flowing from massive spinning particle $i$ to $j (\neq i)$ in the following way <cit.> (see also <cit.>)

Δ_F^μν,ρσ(ℓ) = -iℓ^2 -iϵ 12[K_(ij)^μν,ρσ(ℓ,p_i,p_j)+G_(ij)^μν,ρσ(ℓ,p_i,p_j)],
K^μν,ρσ_(ij)(ℓ,p_i,p_j)    =    𝒞(ℓ,p_i,p_j) [(p_i-ℓ)^μℓ^ν
+(p_i-ℓ)^νℓ^μ] [(p_j+ℓ)^ρℓ^σ+(p_j+ℓ)^σℓ^ρ] ,

G^μν,ρσ_(ij)(ℓ,p_i,p_j)    =   (η^μρη^νσ+η^μση^νρ-η^μνη^ρσ)-K^μν,ρσ_(ij)(ℓ,p_i,p_j) ,
𝒞(ℓ,p_i,p_j)        = (-1){p_i.(p_i-ℓ)-iϵ}
{p_j.(p_j+ℓ)-iϵ} { ℓ.(ℓ-2p_i)-iϵ}  {ℓ.(ℓ+2p_j)-iϵ}

]  .
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 $\widetilde{\Gamma}^{(N+1)}$, we do not carry out any Grammer-Yennie decomposition as those won't contribute to IR divergences. In (<ref>) 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.

This figure is a Feynman diagrammatic representation of the expression in (<ref>). Solid lines represent the massive spinning particles, dashed lines represent the ingoing virtual graviton with momentum $\ell$ and the arrow in the graviton line represents that it is a K-graviton proportional to $\zeta_\mu \ell_\nu +\zeta_\nu \ell_\mu$. The solid blobs in the RHS represent a new kind of vertices and the Feynman rules for the left blob vertex is $-\kappa\lbrace 2 \zeta.(q+\ell) +\zeta^\nu \ell^b \Sigma_{\nu b}\rbrace$ and for the right blob vertex is $\kappa\lbrace 2\zeta.q -\zeta^\nu \ell^b \Sigma_{\nu b}^T \rbrace$.

Ward identities involving K-graviton: From the definition of K-graviton propagator in (<ref>), it is clear that K-graviton is proportional to a pure gauge of structure $\zeta^\mu\ell^\nu +\zeta^\nu \ell^\mu$, with $\zeta=p_i-\ell$ when it flows from $i$-th leg. Let us study the Ward identity for an off-shell un-amputated three particle Green’s function involving two massive fields and one K-graviton. It is diagrammatically represented in Fig.<ref>.

The LHS of the Fig.<ref> takes the following form:
1q^2 +m^2-iϵ1(q+ℓ)^2+m^2-iϵ(ζ^μℓ^ν+ζ^νℓ^μ)[ Ξ(-q)Γ^(3)_μν(q,-q-ℓ,ℓ)Ξ(-q-ℓ)] .
To compute the above expression we used the result of (<ref>) and then simplified using the identities in (<ref>),(<ref>) and derivatives of (<ref>). After all the simplification we get
iκ(q+ℓ)^2+m^2-iϵ{2iζ.ℓΞ(-q)+2iζ.qℓ^ρΞ(-q)q^ρ+iζ.ℓℓ^ρℓ^σ^2 Ξ(-q)q^ρq^σ+2iζ.ℓℓ^ρΞ(-q)q^ρ

+i2ℓ^b ℓ^ρζ^νℓ^μΣ_νb^2 Ξ(-q)q^μq^ρ+i3ζ.q ℓ^ρℓ^σℓ^ν^3 Ξ(-q)q^νq^ρq^σ+2iζ.q Ξ(-q)+iζ^νℓ^bΣ_νbΞ(-q)

+iζ^νℓ^μℓ^b Σ_νbΞ(-q)q^μ+iζ.q ℓ^μℓ^b ^2 Ξ(-q)q^μq^b+𝒪(ℓ^4)}

+iκq^2+m^2-iϵ{-2iζ.q Ξ(-q)+iζ^νℓ^b Ξ(-q)Σ_νb^T } .
Now undoing the $\ell$ expansion for $\Xi(-q-\ell)$ the above expression can be re-written as
1q^2 +m^2-iϵ1(q+ℓ)^2+m^2-iϵ(ζ^μℓ^ν+ζ^νℓ^μ)[ Ξ(-q)Γ^(3)_μν(q,-q-ℓ,ℓ)Ξ(-q-ℓ)]

=   κ{-2 ζ.(q+ℓ) -ζ^νℓ^b Σ_νb}Ξ(-q-ℓ)(q+ℓ)^2+m^2-iϵ+Ξ(-q)q^2+m^2-iϵ κ{2ζ.q -ζ^νℓ^b Σ_νb^T }.

The above identity is diagrammatically represented in Fig.<ref>. 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.<ref>.

In comparison to the QED Ward identity in (<ref>), 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.<ref>, it will result in new Feynman rules for the right blob vertex drawn in the second diagram after the equality of Fig.<ref>. Specifically, the right blob vertex Feynman rule reads $\kappa\lbrace 2\zeta.(q+k) -\zeta^\nu \ell^b \Sigma_{\nu b}^T \rbrace$, which contains an additional term $2\kappa\zeta\cdot k$  relative to the three-point un-amputated Green's function right blob Feynman rule $\kappa\lbrace 2\zeta.q -\zeta^\nu \ell^b \Sigma_{\nu b}^T \rbrace$ in Fig.<ref>. In turn this implies that $\mathcal{A}^{(N+1)}$ contains an extra exponentiation term compared to the IR-exponentiation factor of $\mathcal{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 <cit.>, it is also observed that the right-hand side of Fig.<ref> 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 $\mathcal{A}^{(N+1)}$, which we can systematically account for it as well. However, as mentioned in section-<ref>, we still need to regulate some additional IR divergence in $\mathcal{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
𝒜^(N)≡exp{K_gr} 𝒜^(N)_IR-finite  ,

𝒜^(N+1)reg≡ exp{K_gr} 𝒜^(N+1)_IR-finite ,
K_gr   =     -iκ^2∑_i=1^N  ∑_j=1
j≠i^N ∫d^4ℓ(2π)^4  1ℓ^2-iϵ
1ℓ.(ℓ-2p_i)-iϵ 1ℓ.(ℓ+2p_j)-iϵ

×[2(p_i.p_j)^2-p_i^2p_j^2-ℓ^2(p_i.p_j)+2(p_i.p_j)(p_i.ℓ)-2(p_i.p_j)(p_j.ℓ) ]

is the Eikonal IR-divergent exponentiated factor.
In equations (<ref>) and (<ref>), $\mathcal{A}^{(N)}_{\text{IR-finite}}$ and $\mathcal{A}^{(N+1)}_{\text{IR-finite}}$ represent the infrared finite components of the $N$-particle and $N$-particle-1-graviton amplitudes, respectively.[To determine $\mathcal{A}^{(N+1)}_{\text{IR-finite}}$, we need to use an explicit IR cut-off for the diagrams involving graviton self interaction vertices. The “reg” over the $\equiv$ sign in (<ref>) corresponds to this particular IR regularization scheme. Effectively this IR-regularization scheme removes a factor of $\exp\lbrace K_{phase}\rbrace$ from $\mathcal{A}^{(N+1)}$ with an explicit expression for $K_{phase}$ being
K_phase=  iκ^2 ∑_j=1^N(p_j.k)^2∫_0^R^-1d^4ℓ(2π)^4 1ℓ^2-iϵ1k.ℓ+iϵ1p_j.ℓ-iϵ .
] These components are obtained by removing the exponentiated IR-divergent parts from the original divergent amplitudes defined through the relations (<ref>) and (<ref>). Both $\mathcal{A}^{(N)}_{\text{IR-finite}}$ and $\mathcal{A}^{(N+1)}_{\text{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 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 $\mathcal{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 (<ref>), (<ref>) and (<ref>), here we also split the IR-finite parts of the tree and one-loop amplitudes in the following way
𝒜^(N)_IR-finite,0   ≡    Γ^(N) ,

𝒜^(N+1)_IR-finite,0   ≡    Γ^(N+1) ,

𝒜^(N)_IR-finite,1   ≡   [𝒜^(N)_G,1+𝒜^(N)_K-finite,1+𝒜^(N)_self,1+𝒜^(N)_non-div,1] ,

𝒜^(N+1)_IR-finite,1   ≡   [𝒜^(N+1)_G,1+𝒜^(N+1)_K-finite,1+𝒜^(N+1)_3-graviton-reg,1+𝒜^(N+1)_self,1+𝒜^(N+1)_non-div,1].
Above different components of the IR-finite one-loop amplitudes are defined as

* $\mathcal{A}^{(N)}_{G,1}$ corresponds to the diagram in Fig.<ref> with the dashed line being  a G-graviton propagator.

* $\mathcal{A}^{(N)}_{K-\text{finite},1}$ corresponds to the contribution from the diagram in Fig.<ref> , evaluated with K-graviton propagator representing the dashed line there and then subtracted the contribution $K_{gr}\times \Gamma^{(N)}$ from it at the integrand level.

* $\mathcal{A}^{(N)}_{\text{self,1}}$ represents the set of Feynman diagrams in Fig.<ref> evaluated with full graviton propagator representing the dashed lines.

* $\mathcal{A}^{(N)}_{\text{non-div,1}}$ represents the set of diagrams in Fig.<ref>, evaluated with full graviton propagator representing the dashed lines.

* $\mathcal{A}^{(N+1)}_{G,1}$ corresponds to the set of diagrams in Fig.<ref> 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$.

* $\mathcal{A}^{(N+1)}_{K-\text{finite},1}$ corresponds to the contribution from the diagrams in Fig.<ref>, evaluated with K-graviton propagator representing the dashed line and then subtracted the contribution $K_{gr}\times \Gamma^{(N+1)}$ at the integrand level. The dashed external lines represent the on-shell graviton with momentum $k$.

* $\mathcal{A}^{(N+1)}_\text{3-graviton-reg,1}$ represents the sum of the contributions of Feynman diagrams in Fig.<ref>, when we evaluate them using full graviton propagator and regulate the IR divergence considering detector resolution as the explicit IR cut-off.

* $\mathcal{A}^{(N+1)}_{\text{self,1}}$ represents the set of Feynman diagrams in Fig.<ref>, evaluated with full graviton propagator representing the dashed lines. The dashed external lines represent the on-shell graviton with momentum $k$.

* $\mathcal{A}^{(N+1)}_{\text{non-div,1}}$ represents the set of diagrams in Fig.<ref>, 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 $\omega^{-1}$ and $\ln\omega$ soft factors from the ratio of $\mathcal{A}^{(N+1)}_{\text{IR-finite}}$ and $\mathcal{A}^{(N)}_{\text{IR-finite}}$ when the external graviton energy is small i.e. $\omega<<|p_i^\mu |$.

§.§.§ IR-finite one-loop N-particle amplitude

Let us evaluate the sum of contributions from Fig.<ref> 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.<ref>, where the dashed line represents a full graviton propagator. Then we need to subtract $K_{gr}\times \Gamma^{(N)}$ from the evaluated result. The mathematical expression reads
=    ∑_i,j=1
i> j^N ∫d^4ℓ(2π)^4  1(p_i-ℓ)^2+m_i^2-iϵ 1(p_j+ℓ)^2+m_j^2 -iϵ

×[ ϵ_i^T(-p_i)Γ^(3)_μν (p_i,-p_i+ℓ,-ℓ)Ξ_i(-p_i+ℓ)]Δ_F^μν,ρσ(ℓ)

[ ϵ_j^T(-p_j)Γ^(3)_ρσ (p_j,-p_j-ℓ,ℓ)Ξ_j(-p_j-ℓ)] Γ^(N)_(ij)(p_i-ℓ,p_j+ℓ)

- K_gr×Γ^(N) .
We evaluate the above expression by using the identity derived in (<ref>). The result in small $\ell$ expansion in the integrand turns out to be

=    -iκ^2 ∑_i=1^N∑_j=1
j≠i^Nϵ_i^Tϵ_j^T∫d^4ℓ(2π)^41ℓ^2-iϵ 1ℓ·(ℓ-2p_i)-iϵ 1ℓ·(ℓ+2p_j)-iϵ

[{2(p_i.p_j)^2-p_i^2 p_j^2}(-ℓ^ρp_i^ρΓ^(N)_(ij)(p_i,p_j) +ℓ^ρp_j^ρΓ^(N)_(ij)(p_i,p_j) )

+{2(p_i.p_j)(p_j^ρΣ_iρσ^T ℓ^σ-p_i^ρΣ_jρσ^T ℓ^σ)+p_i^2 p_j^ρΣ_jρσ^T ℓ^σ-p_j^2 p_i^ρΣ_iρσ^T ℓ^σ}Γ^(N)_(ij)(p_i,p_j)

Note that in the limit $\ell^\mu\rightarrow 0$ the integrand of the above expression behaves like $\int \f{d^4\ell}{|\ell|^3}$ at leading order, hence the contribution is IR-finite as promised.

Diagrams in Fig.<ref> 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 (<ref>), the sum of the contribution vanishes i.e.
𝒜^(N)_self,1=0 .
Now let us analyze the diagrams in Fig.<ref> with dashed lines being full graviton propagators are connected to some internal massive virtual lines or massive EFT vertices inside $\Gamma^{(N)}$. In the limit $\ell^\mu\rightarrow 0$ the integrand of the first diagram behaves like $\int \f{d^4\ell}{|\ell|^3}$ at leading order, and the second diagram behaves like $\int \f{d^4\ell}{|\ell|^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
𝒜^(N)_non-div,1 .

Hence the total IR-finite contribution to $N$-particle amplitude follows from the definition (<ref>) becomes
𝒜^(N)_IR-finite,1   =    -iκ^2 ∑_i=1^N∑_j=1
j≠i^Nϵ_i^Tϵ_j^T∫d^4ℓ(2π)^41ℓ^2-iϵ 1ℓ·(ℓ-2p_i)-iϵ 1ℓ·(ℓ+2p_j)-iϵ

[{2(p_i.p_j)^2-p_i^2 p_j^2}(-ℓ^ρp_i^ρΓ^(N)_(ij)(p_i,p_j) +ℓ^ρp_j^ρΓ^(N)_(ij)(p_i,p_j) )

+{2(p_i.p_j)(p_j^ρΣ_iρσ^T ℓ^σ-p_i^ρΣ_jρσ^T ℓ^σ)+p_i^2 p_j^ρΣ_jρσ^T ℓ^σ-p_j^2 p_i^ρΣ_iρσ^T ℓ^σ}Γ^(N)_(ij)(p_i,p_j)

+𝒪(ℓℓ)]  + 𝒜^(N)_non-div,1.

§.§.§ IR-finite one-loop (N+1)-particle amplitude in the soft limit

Here we analyze all the Feynman diagrams contributing to (<ref>) in the soft limit i.e. $\omega\rightarrow 0$. We start by analyzing the diagrams in Fig.<ref> with the dashed lines being full graviton propagators, which evaluates the contribution of $\mathcal{A}^{(N+1)}_{G,1}+
\mathcal{A}^{(N+1)}_{K-\text{finite},1}+K_{gr}\Gamma^{(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 (<ref>). Then finally we subtract the contribution of $K_{gr}\Gamma^{(N+1)}$ from the sum to extract the $\omega^{-1}$ and $\ln\omega$ soft factors from $\mathcal{A}^{(N+1)}_{G,1}+
\mathcal{A}^{(N+1)}_{K-\text{finite},1}$ in the limit $\omega<<|p_i^\mu|$.

The first diagram in Fig.<ref> with full graviton propagator representing the internal dashed line takes the following form
B_I   ≡    ∑_i=1^N∑_j=1
j≠i^Nε^λτ(k)2p_i.k∫d^4ℓ(2π)^4 1(p_i+k-ℓ)^2+m_i^2-iϵ1(p_j+ℓ)^2+m_j^2 -iϵ

[ ϵ_i^T(-p_i)Γ^(3)_λτ (p_i,-p_i-k,k)Ξ_i(-p_i-k)Γ^(3)_ρσ(p_i+k,-p_i-k+ℓ,-ℓ)Ξ_i(-p_i-k+ℓ)]

×Δ_F^ρσ,μν(ℓ) [ϵ_j^T(-p_j)Γ^(3)_μν (p_j,-p_j-ℓ,ℓ)Ξ_j(-p_j-ℓ)] Γ^(N)_(ij)(p_i+k-ℓ,p_j+ℓ) .
After using the identities in (<ref>) and (<ref>) and simplifying the above expression reduces to
B_I   =    iκ^3 ∑_i=1^N∑_j=1
j≠i^Nϵ_i^Tϵ_j^T1p_i.k∫d^4ℓ(2π)^4 1ℓ^2-iϵ1ℓ.(ℓ-2p_i-2k)+2p_i.k -iϵ1ℓ.(ℓ+2p_j)-iϵ

[ (-2(p_i.ε.p_i){2(p_i.p_j)^2-p_i^2p_j^2 +4p_i.p_j p_j.k-2p_i.p_j p_j.ℓ+2p_i.p_j p_i.ℓ}

+4(p_i.ε.p_i) p_i.p_j p_i^μΣ_jμν^Tℓ^ν-2(p_i.ε.p_i) p_i^2 p_j^μΣ_jμν^T ℓ^ν+2{2(p_i.p_j)^2-p_i^2p_j^2}p_iμε^μρΣ_iρν^Tk^ν

-4(p_i.ε.p_i) p_i.p_j p_j^μΣ_iμν^Tℓ^ν+2(p_i.ε.p_i) p_j^2 p_i^μΣ_jμν^T ℓ^ν-2i p_i.k p_i.p_j p_iμε^μν𝒦_i(-p_i)p_i^νΞ_i(-p_i)p_i^σp_j^σ

+ ip_i.k p_j^2 p_iμε^μν𝒦_i(-p_i)p_i^νΞ_i(-p_i)p_i^σp_i^σ)Γ^(N)_(ij)(p_i,p_j) -2 (p_i.ε.p_i){2(p_i.p_j)^2-p_i^2p_j^2 }

×((k-ℓ)^μp_i^μΓ^(N)_(ij)(p_i,p_j) +ℓ^μp_j^μΓ^(N)_(ij)(p_i,p_j)) + 𝒪(ℓℓ,ℓk, kk)] .
The second diagram in Fig.<ref> with full graviton propagator representing the dashed internal line takes the following form
B_II   ≡    ∑_i=1^N∑_j=1
j≠i^N∫d^4ℓ(2π)^4 1(p_i+k-ℓ)^2+m_i^2-iϵ1(p_j+ℓ)^2+m_j^2 -iϵ Δ_F^ρσ,λτ(ℓ)

×[ ϵ_i^T(-p_i)ε^μν(k)Γ^(4)_μν, ρσ (p_i,-p_i-k+ℓ,k,-ℓ)Ξ_i(-p_i-k+ℓ)]

×  [ϵ_j^T(-p_j)Γ^(3)_λτ (p_j,-p_j-ℓ,ℓ)Ξ_j(-p_j-ℓ)] Γ^(N)_(ij)(p_i+k-ℓ,p_j+ℓ) .
After using the identities in (<ref>) and (<ref>) and simplifying the above expression reduces to
B_II    =    iκ^3 ∑_i=1^N∑_j=1
j≠i^Nϵ_i^Tϵ_j^T∫d^4ℓ(2π)^4 1ℓ^2-iϵ1ℓ.(ℓ-2p_i-2k)+2p_i.k -iϵ1ℓ.(ℓ+2p_j)-iϵ

[ {-4(p_i.ε.p_i)p_j^2+16 (p_i.ε.p_j)(p_i.p_j)+2ip_i.p_j p_iμε^μν𝒦_i(-p_i)p_i^νp_j^σΞ_i(-p_i)p_i^σ

+2ip_i.p_j p_j^σ𝒦_i(-p_i)p_i^σ p_iμε^μνΞ_i(-p_i)p_i^ν -ip_j^2 p_iμε^μν𝒦_i(-p_i)p_i^νp_i^σΞ_i(-p_i)p_i^σ

-ip_j^2 p_i^σ𝒦_i(-p_i)p_i^σp_iμε^μνΞ_i(-p_i)p_i^ν}Γ^(N)_(ij)(p_i,p_j) +𝒪(ℓ, k)] .
The third diagram in Fig.<ref> with full graviton propagator representing the dashed internal line takes the following form
B_III   ≡    ∑_i=1^N∑_j=1
j≠i^N∫d^4ℓ(2π)^4  1(p_i-ℓ)^2+m_i^2-iϵ1(p_i+k-ℓ)^2+m_i^2-iϵ1(p_j+ℓ)^2+m_j^2 -iϵ

[ ϵ_i^T(-p_i)Γ^(3)_μν (p_i,-p_i+ℓ,-ℓ)Ξ_i(-p_i+ℓ)ε^ρσ(k)Γ^(3)_ρσ(p_i-ℓ,-p_i-k+ℓ,k)Ξ_i(-p_i-k+ℓ)]

×Δ_F^μν,λτ(ℓ)  [ϵ_j^T(-p_j)Γ^(3)_λτ (p_j,-p_j-ℓ,ℓ)Ξ_j(-p_j-ℓ)] Γ^(N)_(ij)(p_i+k-ℓ,p_j+ℓ) .
After using the identities in (<ref>) and (<ref>) and simplifying the above expression reduces to
B_III   =    iκ^3 ∑_i=1^N∑_j=1
j≠i^Nϵ_i^Tϵ_j^T∫d^4ℓ(2π)^4 1ℓ^2-iϵ1ℓ.(ℓ-2p_i) -iϵ1ℓ.(ℓ-2p_i-2k)+2p_i.k -iϵ

×1ℓ.(ℓ+2p_j)-iϵ[ (-4(p_i.ε.p_i){2(p_i.p_j)^2 -p_i^2 p_j^2 }-8(p_i.ε.p_i) p_i.p_j p_i.ℓ

+8 (p_i.ε.ℓ) {2(p_i.p_j)^2 -p_i^2 p_j^2 }+8 (p_i.ε.p_i) p_i.p_j p_j.ℓ+4 {2(p_i.p_j)^2 -p_i^2 p_j^2 }p_iρε^ρσΣ_iσμ^Tk^μ

-8p_i.p_j (p_i.ε.p_i) p_j^μΣ_iμν^T ℓ^ν+4p_j^2 (p_i.ε.p_i) p_i^μΣ_iμν^T ℓ^ν+8(p_i.ε.p_i) p_i.p_j p_i^μΣ_jμν^T ℓ^ν

-4p_i^2 (p_i.ε.p_i) p_j^μΣ_jμν^T ℓ^ν+4i p_i.ℓp_i.p_j p_j^μ𝒦_i(-p_i)p_i^μp_iρε^ρσΞ_i(-p_i)p_i^σ

-2i p_i.ℓp_j^2 p_i^μ𝒦_i(-p_i)p_i^μp_iρε^ρσΞ_i(-p_i)p_i^σ)Γ^(N)_(ij)(p_i,p_j) -4 (p_i.ε.p_i){2(p_i.p_j)^2-p_i^2p_j^2 }

×((k-ℓ)^μp_i^μΓ^(N)_(ij)(p_i,p_j) +ℓ^μp_j^μΓ^(N)_(ij)(p_i,p_j)) + 𝒪(ℓℓ,ℓk, kk)] .
The fourth diagram in Fig.<ref> with full graviton propagator representing the dashed internal line takes the following form
B_IV   ≡    ∑_k=1^N 12p_k·kϵ_k^T(-p_k)ε^λτ(k)Γ^(3)_λτ (p_k,-p_k-k,k)Ξ_k(-p_k-k)

×12 ∑_i=1
j≠i,k^N ∫d^4ℓ(2π)^4  1(p_i-ℓ)^2+m_i^2-iϵ 1(p_j+ℓ)^2+m_j^2 -iϵ

×[ ϵ_i^T(-p_i)Γ^(3)_μν (p_i,-p_i+ℓ,-ℓ)Ξ_i(-p_i+ℓ)]Δ_F^μν,ρσ(ℓ)

×[ ϵ_j^T(-p_j)Γ^(3)_ρσ (p_j,-p_j-ℓ,ℓ)Ξ_j(-p_j-ℓ)] Γ^(N)_(ijk)(p_i-ℓ,p_j+ℓ,p_k+k) .
The fifth diagram in Fig.<ref> with full graviton propagator representing the dashed internal line takes the following form
B_V   ≡     12∑_i=1^N∑_j=1
j≠i^N ∫d^4ℓ(2π)^4  1(p_i-ℓ)^2+m_i^2-iϵ 1(p_j+ℓ)^2+m_j^2 -iϵ

×[ ϵ_i^T(-p_i)Γ^(3)_μν (p_i,-p_i+ℓ,-ℓ)Ξ_i(-p_i+ℓ)]Δ_F^μν,ρσ(ℓ) ε^λτ(k)

×[ ϵ_j^T(-p_j)Γ^(3)_ρσ (p_j,-p_j-ℓ,ℓ)Ξ_j(-p_j-ℓ)] Γ^(N+1)_(ij), λτ( p_i-ℓ,p_j+ℓ,k) .
In both the expressions of $B_{IV}$ and $B_{V}$ the integrands can be simplified using (<ref>) and the result of the following common loop integral reads
∫d^4ℓ(2π)^4  1(p_i-ℓ)^2+m_i^2-iϵ 1(p_j+ℓ)^2+m_j^2 -iϵ[ ϵ_i^T(-p_i)Γ^(3)_μν (p_i,-p_i+ℓ,-ℓ)Ξ_i(-p_i+ℓ)]

×Δ_F^μν,ρσ(ℓ) ×[ ϵ_j^T(-p_j)Γ^(3)_ρσ (p_j,-p_j-ℓ,ℓ)Ξ_j(-p_j-ℓ)]

=   -2iκ^2 ϵ_i^Tϵ_j^T∫d^4ℓ(2π)^41ℓ^2-iϵ 1ℓ·(ℓ-2p_i)-iϵ 1ℓ·(ℓ+2p_j)-iϵ[  2(p_i.p_j)^2-p_i^2p_j^2

+2(p_i.p_j)(p_i.ℓ)-2(p_i.p_j)(p_j.ℓ)  +2(p_i.p_j)(p_j^ρΣ_iρσ^T ℓ^σ-p_i^ρΣ_jρσ^T ℓ^σ)+p_i^2 p_j^ρΣ_jρσ^T ℓ^σ

-p_j^2 p_i^ρΣ_iρσ^T ℓ^σ +𝒪(ℓℓ)].
Hence the total IR-finite contribution from the sets of diagrams in Fig.<ref> reads

=    B_I+B_II+B_III+B_IV+B_V-K_gr×Γ^(N+1) ,
where the $(N+1)$ particle tree level amplitude is given by
Γ^(N+1)    =    ∑_k=1^N 12p_k·kϵ_k^T(-p_k)ε^λτ(k)Γ^(3)_λτ (p_k,-p_k-k,k)Ξ_k(-p_k-k) Γ^(N)_(k)(p_k+k)

+ ε^λτ(k) Γ^(N+1)_λτ(k) .
In the soft limit the above tree level $(N+1)$ particle amplitude provides the tree level soft graviton theorem with soft factor given in (<ref>).

Substituting the results of (<ref>), (<ref>), (<ref>), (<ref>), (<ref>) and (<ref>) in (<ref>), it is easy to see that the final expression is IR-finite in the limit $\ell^\mu\rightarrow 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}\times \Gamma^{(N+1)}$, and $B_{II}$ is IR finite at finite value of $k$.

Now let us analyze the expression (<ref>) in the soft limit i.e. $\omega<<|p_i^\mu|$. In this limit a part of the final expression (<ref>) contributes to $\omega^{-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 (<ref>) we divide the virtual momenta integration range into three regions: $|\ell^\mu|\in [0,\omega], [\omega ,|p_i^\mu|]$ and $[|p_i^\mu|,\infty)$. It turns out that the integrand starts contributing at order $\omega^0$ in the region of integration $|\ell^\mu|\in [0,\omega]$ and $[|p_i^\mu|,\infty)$[In the integration region $|\ell^\mu|\in [|p_i^\mu|,\infty)$ the integrand is UV divergent, so needs to use UV regulator and add appropriate counter terms to extract finite result. ] in the soft limit. Finally in the region of integration $|\ell^\mu|\in [\omega ,|p_i^\mu|]$ the integrand contributes at order $\ln\omega$ which is dominant compare to order $\omega^0$ in the soft expansion.  Hence the order $\omega^{-1}$ and $\ln\omega$ contribution from (<ref>) turns out to be

=     κ∑_k=1^N p_k. ε. p_kp_k. k ×(𝒜^(N)_G,1+𝒜^(N)_K-finite,1)

-iκ^34 ∑_i=1^N∑_j=1
j≠i^Nϵ_i^Tϵ_j^T∫_reg d^4ℓ(2π)^4 1ℓ^2-iϵ 1ℓ.p_i +iϵ1ℓ.p_j-iϵ

×[-8(p_i.ε.p_i) p_i.p_j p_j.kp_i.k  +2ℓ.kp_i.k1ℓ.p_i +iϵ (p_i.ε.p_i){2(p_i.p_j)^2-p_i^2p_j^2 }

- 4 (p_i.ε.ℓ)ℓ.p_i+iϵ   {2(p_i.p_j)^2 -p_i^2 p_j^2 }-4(p_i.ε.p_i)p_j^2+16 (p_i.ε.p_j)(p_i.p_j)] Γ^(N)_(ij)(p_i,p_j)

+𝒪(ω^0) .
Above the expression of $\mathcal{A}^{(N)}_{G,1}+\mathcal{A}^{(N)}_{K-\text{finite},1}$ is given in (<ref>). The “$reg$” in the subscript of loop momentum integration above refers to the integration range $|\ell^\mu|\in [\omega ,|p_i^\mu|],|p_j^\mu|]$. It is intriguing to observe that the components reliant on theory, such as $\mathcal{K}_i$, $\Xi_i$, and $\frac{\partial\Gamma^{(N)}}{\partial 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\omega$ order soft factor contributing integrand. This theory and spin independence feature is not true for the order $\omega^0$ contribution from (<ref>), which we ignored here.

Diagrams containing 3-graviton vertex contributing to $\mathcal{A}^{(N+1)}_\text{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.

Let us proceed to calculate the contribution of $\mathcal{A}^{(N+1)}_\text{3-graviton-reg,1}$. This term represents the sum of contributions from the Feynman diagrams shown in Fig.<ref>.  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.<ref> reads
≡    ∑_i,j=1
i> j^N ∫_R^-1^∞ d^4ℓ(2π)^4  1(p_i-ℓ)^2+m_i^2-iϵ 1(p_j+ℓ+k)^2+m_j^2 -iϵ

×[ ϵ_i^T(-p_i)Γ^(3)_ab (p_i,-p_i+ℓ,-ℓ)Ξ_i(-p_i+ℓ)]

×Δ_F^ab,ρσ(ℓ)ε^μν(k)V^(3)_μν, ρσ,τκ(k,ℓ, -ℓ-k) Δ_F^τκ,cd(ℓ)

[ ϵ_j^T(-p_j)Γ^(3)_cd (p_j,-p_j-ℓ-k,ℓ+k)Ξ_j(-p_j-ℓ-k)] Γ^(N)_(ij)(p_i-ℓ,p_j+ℓ+k)

Now using the identity in (<ref>) and substituting the 3-graviton vertex from (<ref>) with on-shell and transverse-traceless condition for the external graviton the above expression simplifies to
=   -iκ^38 ∑_i,j=1
i≠j^N ϵ_i^T ϵ_j^T∫_R^-1^∞d^4ℓ(2π)^4  1(p_i-ℓ)^2+m_i^2-iϵ 1(p_j+ℓ+k)^2+m_j^2 -iϵ1ℓ^2-iϵ

×1(ℓ+k)^2 -iϵ[4p_i^ρp_i^σ-2p_i^2η^ρσ -2p_i^ρℓ^σ-2p_i^σℓ^ρ+2p_i^ρΣ^Tσa_iℓ_a+2p_i^σΣ^Tρa_iℓ_a

-2η^ρσp_i^aΣ^T_iacℓ^c +𝒪(ℓℓ)  ]×ε^μν(k)[  -(k.ℓ+ℓ^2)(

+ 4(k.ℓ+ℓ^2)η_ρνη_στη_κμ- 2  (k_τℓ_κη_μρη_νσ-ℓ_μℓ_νη_ρτη_σκ-(k+ℓ)_ρ k_ση_μτη_νκ)

-2  {(ℓ_μℓ_τ-k_τℓ_μ)η_νρη_σκ+(2ℓ_μk_ρ+ℓ_μℓ_ρ)η_τνη_κσ+(k_ρk_τ +ℓ_ρk_τ-ℓ_τk_ρ)η_μση_νκ}

+𝒪(k^4, k^3ℓ, k^2ℓ^2, kℓ^3,ℓ^4)]×[4p_j^τp_i^κ-2p_j^2η^τκ +2p_j^τ(ℓ+k)^κ+2p_j^κ(ℓ+k)^τ

-2p_j^τΣ^Tκb_j(ℓ+k)_b-2p_j^κΣ^Tτb_j(ℓ+k)_b +2η^τκp_j^bΣ^T_jbd(ℓ+k)^d  +𝒪(ℓℓ,ℓk, kk)]

×  Γ^(N)_(ij)(p_i-ℓ,p_j+ℓ+k)  .
The above expression contributes at order $\ln\omega$ in three regions of integration $|\ell^\mu|\in [R^{-1},\omega]$, $|(\ell+k)^\mu|\in [R^{-1},\omega]$, and  “$reg$”$\equiv |\ell^\mu|\in [\omega, |p_i^\mu|]$ in the soft limit. The contribution from the regions $|\ell^\mu|\in [R^{-1},\omega]$ and $|(\ell+k)^\mu|\in [R^{-1},\omega]$ turn out to be same as $C_{I}$ is symmetric under the simultaneous exchange of $\ell\leftrightarrow (\ell+k)$ and $i\leftrightarrow j$. Hence we only evaluate the contribution in region $|\ell^\mu|\in [R^{-1},\omega]$ and multiply by a factor of $2$ to incorporate the contribution of other region.

In the region $|\ell^\mu|\in [R^{-1},\omega]$ we approximate the propagator denominators at leading order as
1(p_i-ℓ)^2+m_i^2-iϵ≃-12p_i.ℓ+iϵ  ,   1(p_j+ℓ+k)^2+m_j^2 -iϵ≃12p_j.k-iϵ  ,

1(ℓ+k)^2 -iϵ ≃12ℓ.k-iϵ  .
Then only when we choose the order $\mathcal{O}(kk)$ terms from the $V^{(3)}$ vertex, after the loop momentum integration we can have $\mathcal{O}(\ln\omega)$ contribution. On the other hand in the region “$reg$”$\equiv |\ell^\mu|\in [\omega, |p_i^\mu|]$ we approximate the propagator denominators at leading order as
1(p_i-ℓ)^2+m_i^2-iϵ≃-12p_i.ℓ+iϵ  ,   1(p_j+ℓ+k)^2+m_j^2 -iϵ≃12p_j.ℓ-iϵ  ,

1(ℓ+k)^2 -iϵ ≃1ℓ^2-iϵ  .
Then only when we choose the order $\mathcal{O}(\ell\ell)$ terms from the $V^{(3)}$ vertex, after the loop momentum integration we can have $\mathcal{O}(\ln\omega)$ contribution.

Hence the full contribution at order $\ln\omega$ from all the three regions of integration becomes
=   iκ^34 ∑_i,j=1
i≠j^Nϵ_i^T ϵ_j^T1p_j.k∫_R^-1^ωd^4ℓ(2π)^4  1ℓ^2-iϵ1p_i.ℓ+iϵ1ℓ.k-iϵ[ 4(p_i.k)^2 (p_j.ε.p_j)

-4(p_i.k)(p_j.k)(p_i.ε.p_j)] Γ^(N)_(ij)(p_i,p_j)

+iκ^32 ∑_i,j=1
i≠j^N ϵ_i^T ϵ_j^T∫_reg d^4ℓ(2π)^4  1{ℓ^2-iϵ}^21p_i.ℓ+iϵ1p_j.ℓ-iϵ[ -ℓ^2 (p_i.ε.p_i)p_j^2

-ℓ^2 (p_j.ε.p_j)p_i^2 +4ℓ^2 (p_i.p_j)(p_i.ε.p_j) +2(p_i.p_j)^2 (ℓ.ε.ℓ)-2(p_i.ε.ℓ)(p_i.p_j)(p_j.ℓ)

+2p_j^2 (p_i.ε.ℓ)(p_i.ℓ)-p_i^2 p_j^2 (ℓ.ε.ℓ) -2(p_j.ε.ℓ)(p_i.p_j)(p_i.ℓ)+2p_i^2 (p_j.ε.ℓ)(p_j.ℓ)]

×  Γ^(N)_(ij)(p_i,p_j)  +  𝒪(ω^0) ,
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 $\mathcal{K}_i, \Xi_i$ or non-minimal couplings. The first integrand above can be evaluated using the result of the integral (derived in <cit.>)
∫_R^-1^ωd^4ℓ(2π)^4  1ℓ^2-iϵ1p_i.ℓ+iϵ1ℓ.k-iϵ

=    -14π1p_i.kln(ωR)[δ_η_i,-1 -i2π ln(p_i^2(p_i.𝐧)^2)]+𝒪(ω^-1) .
and momentum conservation relation
j≠i^N p_j^μ=-p_i^μ .
Above $\eta_i$ convention is the same as described below (<ref>).
We won't evaluate the second integrand in (<ref>) explicitly at this moment, but will simplify some of the terms which contains $(p_i.\ell)$ or $(p_j.\ell)$ in the numerator using the momentum conservation relation in (<ref>). Also the term containing $(\ell.\varepsilon.\ell)$ can be simplified by the following using integration by parts relation
∫d^4ℓ(2π)^4ℓ^μℓ^ν(ℓ^2 -iϵ)^21p_i.ℓ+iϵ1p_j.ℓ-iϵ

=   12∫d^4ℓ(2π)^41ℓ^2-iϵ1p_i.ℓ+iϵ1p_j.ℓ-iϵ(η^μν-p_i^μℓ^νp_i.ℓ+iϵ-p_j^μℓ^νp_j.ℓ-iϵ) .
Following all the steps outlined above and using exchange symmetry $(p_i,\ell)\leftrightarrow (p_j,-\ell)$ in the second integrand of (<ref>), the simplified expression of $C_I$ becomes
=   -iκ^34π lnω ∑_j=1^Nϵ_j^Tp_j.ε.p_jp_j.k×∑_i=1^N ϵ_i^T  (p_i.k)[ δ_η_i,-1 -i2π ln(p_i^2(p_i.𝐧)^2)]  Γ^(N)_(ij)(p_i,p_j)

+iκ^32 ∑_i,j=1
i≠j^N ϵ_i^T ϵ_j^T∫_reg d^4ℓ(2π)^4  1ℓ^2-iϵ1p_i.ℓ+iϵ1p_j.ℓ-iϵ[ - (p_i.ε.p_i)p_j^2 - (p_j.ε.p_j)p_i^2

+4 (p_i.p_j)(p_i.ε.p_j) -1p_i.ℓ+iϵ{2(p_i.p_j)^2 -p_i^2 p_j^2 }(p_i.ε.ℓ)]Γ^(N)_(ij)(p_i,p_j)   + 𝒪(ω^0) .
Note that in the first line above $j=i$ sum is included while it was not present in the expression (<ref>). The inclusion of this term originates from the second term in the numerator of the first integral in (<ref>) after using the momentum conservation relation (<ref>).

For completeness here we also briefly analyze the second diagram in Fig.<ref>  which reads
C_II   ≡   ∑_i=1^N 12p_i.kϵ_i^T∫_R^-1^∞ d^4ℓ(2π)^41(p_i-ℓ)^2 +m_i^2-iϵ   Γ^(3)_ab(p_i,-p_i+ℓ,-ℓ)Ξ_i(-p_i+ℓ)

×Γ^(3)_cd(p_i-ℓ, -p_i-k,ℓ+k)Ξ_i(-p_i-k)Δ_F^ab,ρσ(ℓ)ε^μνV^(3)_μν,ρσ,τκ(k,ℓ,-k-ℓ)

×Δ_F^τκ,cd(ℓ+k)Γ^(N)_(i)(p_i+k) .
In the integration region  $|\ell^\mu|\in [R^{-1},\omega]$ or $|(\ell+k)^\mu|\in [R^{-1},\omega]$ the numerator of the potentially contributing $\ln\omega$ 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\omega$ becomes
C_II   =    -iκ^32∑_i=1^Nϵ_i^T∫_reg d^4ℓ(2π)^4 1(p_i.ℓ+iϵ)^21ℓ^2 -iϵ [-2p_i^2 (p_i.ε.p_i)+(p_i^2)^2p_i.ℓ+iϵ (p_i.ε.ℓ)]

×Γ^(N)_(i)(p_i) + 𝒪(ω^0) .
The above expression of $C_{II}$ can be simplified using the following identity
∫_reg d^4ℓ(2π)^4 ℓ^μ(p_i.ℓ+iϵ)^31ℓ^2 -iϵ=1p_i^2 ∫_reg d^4ℓ(2π)^4 1(p_i.ℓ+iϵ)^21ℓ^2 -iϵ ,
and the simplified expression reads
C_II   =    iκ^32∑_i=1^Nϵ_i^T∫_reg d^4ℓ(2π)^4 1(p_i.ℓ+iϵ)^21ℓ^2 -iϵ p_i^2 (p_i.ε.p_i)×Γ^(N)_(i)(p_i) +𝒪(ω^0) .

The third, fourth and fifth diagrams in Fig.<ref>  do not contribute to order $\ln\omega$ from any integration region. Hence the total order $\ln\omega$ contribution to $\mathcal{A}^{(N+1)}_\text{3-graviton-reg,1}$ after summing over (<ref>) and (<ref>) turns out to be

=   -iκ^34π lnω ∑_j=1^Nϵ_j^Tp_j.ε.p_jp_j.k×∑_i=1^N ϵ_i^T  (p_i.k)[ δ_η_i,-1 -i2π ln(p_i^2(p_i.𝐧)^2)]  Γ^(N)_(ij)(p_i,p_j)

+iκ^32 ∑_i,j=1
i≠j^N ϵ_i^T ϵ_j^T∫_reg d^4ℓ(2π)^4  1ℓ^2-iϵ1p_i.ℓ+iϵ1p_j.ℓ-iϵ[ - (p_i.ε.p_i)p_j^2 - (p_j.ε.p_j)p_i^2

+4 (p_i.p_j)(p_i.ε.p_j) -1p_i.ℓ+iϵ{2(p_i.p_j)^2 -p_i^2 p_j^2 }(p_i.ε.ℓ)]Γ^(N)_(ij)(p_i,p_j)

+iκ^32∑_i=1^Nϵ_i^T∫_reg d^4ℓ(2π)^4 1(p_i.ℓ+iϵ)^21ℓ^2 -iϵ p_i^2 (p_i.ε.p_i)×Γ^(N)_(i)(p_i)  + 𝒪(ω^0) .

Set of 1-loop diagrams contributing to $\mathcal{A}^{(N+1)}_{\text{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.

Let us now analyze the set of Feynman diagrams in Fig.<ref> which contributes to $\mathcal{A}^{(N+1)}_{\text{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\omega$ 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  (<ref>). 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.<ref>. 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\omega$ in the soft limit
𝒜^(N+1)_self,1   =    -iκ^32 ∑_i=1^Nϵ_i^T∫_reg d^4ℓ(2π)^4 1(p_i.ℓ+iϵ)^21ℓ^2 -iϵ p_i^2 (p_i.ε.p_i)×Γ^(N)_(i)(p_i)

+𝒪(ω^0) .

Before even evaluating $\mathcal{A}^{(N+1)}_{\text{non-div,1}}$ from the set of diagrams in Fig.<ref>, it becomes evident that these diagrams are IR-finite for finite $k$ when $\ell^\mu\rightarrow 0$.  The sum of first, fourth and sixth diagrams in Fig.<ref> contribute to leading soft graviton factor at order $\omega^{-1}$ multiplying $\mathcal{A}^{(N)}_{\text{non-div,1}}$. On the other hand in the integration region $\omega<<|\ell^\mu|<<|p_i^\mu |$, though individually the first and third diagrams in Fig.<ref> have the potential to contribute at order $\ln\omega$ as those behave like $\int_\omega\f{d^4\ell}{|\ell|^4}$, when we sum them up, in the numerator we get an extra factors of $\ell$ and/or  $k$. Hence they can only contribute from order $\omega^0$ or $\omega\ln\omega$. Hence, summing over all the contribution we get
𝒜^(N+1)_non-div,1=κ∑_i=1^Np_i.ε.p_ip_i.k 𝒜^(N)_non-div,1 +𝒪(ω^0) .

Summing over the contribution of (<ref>), (<ref>), (<ref>) and (<ref>) we get
𝒜^(N+1)_IR-finite,1   =   κ∑_i=1^Np_i.ε.p_ip_i.k ×𝒜^(N)_IR-finite,1

-iκ^34π lnω ∑_j=1^Np_j.ε.p_jp_j.k×∑_i=1^N   (p_i.k)[ δ_η_i,-1 -i2π ln(p_i^2(p_i.𝐧)^2)] ×Γ^(N)

-iκ^32 ∑_i=1^N∑_j=1
j≠i^N∫_reg d^4ℓ(2π)^4 1ℓ^2-iϵ 1ℓ.p_i +iϵ1ℓ.p_j-iϵ

×[-4(p_i.ε.p_i) p_i.p_j p_j.kp_i.k  +ℓ.kp_i.k1ℓ.p_i +iϵ (p_i.ε.p_i){2(p_i.p_j)^2-p_i^2p_j^2 }

-  (p_i.ε.ℓ)ℓ.p_i+iϵ   {2(p_i.p_j)^2 -p_i^2 p_j^2 }+4 (p_i.ε.p_j)(p_i.p_j)]×Γ^(N)  + 𝒪(ω^0) ,
where the expression of $\mathcal{A}^{(N)}_{\text{IR-finite,1}}$ is given in (<ref>).

Soft graviton theorem result at one-loop:
The expression for one-loop amplitude in the soft expansion derived in (<ref>) can be re-written in the following compact way
𝒜^(N+1)_IR-finite,1   =    κ∑_i=1^Np_i.ε.p_ip_i.k ×𝒜^(N)_IR-finite,1

+  K_phase^reg×κ∑_i=1^Np_i.ε.p_ip_i.k ×𝒜^(N)_IR-finite,0

+κ∑_i=1^Nε_μνp_i^μk_ρp_i·k {p_i^νK_gr^regp_iρ-p_i^ρK_gr^regp_iν}×𝒜^(N)_IR-finite,0   + 𝒪(ω^0) ,
where the expressions of $K_{phase}^{reg}$ and $K_{gr}^{reg}$ are given in (<ref>) and (<ref>) respectively. This result agrees with the loop corrected subleading soft graviton theorem, originally derived in <cit.> for minimally coupled scalar-gravity with scalar contact interaction. Here the re-derivation 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\omega$ 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\omega$ soft graviton factor in (<ref>) undergoes correction due to electromagnetic interactions. The correction term has been derived in <cit.> in a minimally coupled charged scalar theory. A straightforward generalization of our derivation, combined with the covariantization prescription outlined in <cit.>, 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 $\omega\ln\omega$ soft graviton factor as provided in (<ref>), once we understand how to account for the additional subtleties highlighted in the discussion below (<ref>). Furthermore, through an analysis of two-loop amplitudes in the soft limit, it is possible to derive the order $\omega(\ln\omega)^2$ soft graviton theorem as conjectured in (<ref>), which is also expected to be universal. We intend to pursue these analyses in the future.

Using the Feynman rules derived in section-<ref>, in combination with the identities in appendix-<ref>, it is possible to derive the tree-level simultaneous and consecutive double soft-graviton 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\omega$ 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.[BS acknowledges Waël Aoun for pursuing his master thesis project on this topic using the IR-finite S-matrix formalism proposed in <cit.> and making substantial progress.] 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-<ref>, with some of these theorems conjectured based on classical analysis.

In <cit.>, 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 [The explicit expression is given in equation (2.5) of <cit.>, 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.]
ε^μν(k)e_μν(ω, Rn)   =    -i 2GR exp{iω( R+2GlnR∑_j=1
η_j=-1^N p_j·𝐧)}×𝕊^gr_classical ,
where $\mathbb{S}^{gr}_{\text{classical}}$ represents the classical limit of quantum soft factor for single soft graviton emission with momentum $k^{\mu}=-\omega \mathbf{n}^\mu$. In the above expression $\widetilde{e}_{\mu\nu}$ is defined by
e_μν(ω, x⃗)   ≡    ∫_-∞^∞dt e^iωt (h_μν(t,x⃗)-12η_μνh^ρ_ρ(t,x⃗)) .