Even-point Multi-loop Unitarity and its Applications: Exponentiation, Anomalies and Evanescence

We identify novel structure in newly computed multi-loop amplitudes and quantum actions for even-point effective field theories, including both the nonlinear sigma model (NLSM) and double-copy gauge theories such as Born-Infeld and its supersymmetric generalizations. We exploit special properties of all even-point theories towards efficient unitarity based amplitude construction. In doing so, we find evidence that the leading IR divergence of NLSM amplitudes exponentiates when the symmetry group is $\mathbb{CP}^1\cong SU(2)/U(1)$. We then systematically compute the two-loop anomalous behavior of Born-Infeld, and find that the counterterms needed to restore $U(1)$ invariant behavior at loop-level can be constructed via a symmetric-structure double-copy. We also demonstrate that the divergent part of the one-minus $(-+++)$ two-loop anomaly vanishes upon introducing an evanescent operator. In addition to these pure photon counterterms, we verify through explicit calculation that the anomalous matrix elements that violate $U(1)$ duality invariance can be alternatively cancelled by summing over internal $\mathcal{N}=4$ DBIVA superfields. Finally we find that $\mathcal{N}=4$ Dirac-Born-Infeld-Volkov-Akulov (DBIVA) amplitudes permit double-copy construction through two-loop order by reproducing our unitarity based result with a double copy between color-dual $\mathcal{N}=4$ super-Yang-Mills and our two-loop NLSM amplitudes. This result supports the possibility of color-dual representations for NLSM beyond one-loop. We conclude with an overview of how $D$-dimensional four-photon counterterms can be constructed in generality with the symmetric-structure double-copy, and outline a convenient way of counting evanescent operators using Hilbert series as generating functions.


Introduction
Recent decades have seen significant advances in our ability to compute scattering amplitudes to higher orders in perturbation theory. At the epicenter of this explosion in the generation of sharp S-matrix data is the unitarity method [1][2][3], which bypasses the standard Feynman diagram approach by constructing higher order amplitudes directly from lower order on-shell information. At one-loop, the unitarity method allows one to extract amplitudes by directly computing a series of unitarity cuts [4]. On-shell methods have furthermore unveiled novel amplitude-level structure, like the duality between color and kinematics [5] and associated double-copy construction [6], which are frequently obscured by the traditional Lagrangian description of the theories. While much of the research in perturbative calculations thus far has focused on gauge theory and gravity, recent literature has investigated to what extent the S-matrix of effective field theory is constrained by on-shell data [7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25].
In this paper we make the quantum leap to the multi-loop sector of effective field theory. Namely we study even-point (EP) effective field theory (EFT) at the multi-loop level using a combination of generalized unitarity [1][2][3] and the double-copy construction [5,6]. The EFTs that we consider are the nonlinear sigma model (NLSM) and a related family of gauge theories involving Born-Infeld theory and supersymmetric generalizations known as Dirac-Born-Infeld-Volkov-Akulov (DBIVA) theories. These families are known at tree-level to be double-copies between NLSM and Yang-Mills theories with varying amounts of supersymmetry. While we focus on these specific examples, the methods we develop can be used for perturbative calculations in any even-point effective field theory. The motivation for this work is three-fold.
1. First, perturbative calculations in NLSM and DBIVA effective field theories are in many ways significantly simpler than the typical multi-loop calculations of more phenomenological theories, like quantum-chromodynamics (QCD). As we will show, the multi-loop amplitudes for these theories allow us to recycle many of the D-dimensional integration tools that are so powerful at one-loop order. To approach this problem, we expand on the now standard procedure of Forde [4], that computes one-loop amplitudes directly from 4D unitarity cuts. In section 3 we'll show that there are only two basis integrals needed at two-loop four-point for even-point effective field theory amplitudes. In the spirit of one-loop unitarity, the integral coefficients can be directly extracted from Even-point Multi-loop Unitarity (EMU) and a tensor reduction algorithm. Our method of EMU is a D-dimensional approach to integrand construction. This allows us to compute a large catalog of D-dimensional twoloop amplitudes, which provides fertile ground for cultivating insights about the multi-loop structure of effective field theories. The perturbative depth of our calculations sheds light on the exponential structure of IR divergences for CP 1 NLSM amplitudes through two-loops, the anomalous behavior U (1) duality invariance through two-loop order in DBIVA theories, as well as the emergence of evanescent operators relevant for anomaly cancellation in pure Born-Infeld theory. We believe these surprisingly rich physical structures could serve as a theoretical laboratory for future studies of the interplay between effective field theory operators and perturbative calculations in quantum field theory.
2. The second motivation is that very little is known about the duality between color and kinematics at multi-loop level for generic models in the web of theories [26]. For a comprehensive review of color-dual representations, see Refs. [26][27][28] and references therein. While there has been tremendous success of applying the color-kinematics duality to N = 4 super-Yang-Mills, for which color-dual integrands are known through four-point four-loop [29][30][31], five-point through three-loops, and to seven-point at one-loop [32][33][34], there are many obstructions for generic gauge/gravity theories. Presently, the state-of-the-art for nonlinear sigma model (NLSM) color-dual numerators come from the XYZ model of Cheung and Shen [35], while D-dimensional integrands for pure Yang-Mills have only been identified through five-point one-loop [36]. At present, there are no known D-dimensional representations for either of these theories at two-loop that globally manifest the duality between color and kinematics, despite attempts in the literature [37][38][39][40]. Furthermore, similar bottlenecks exists for less-than-maximal N < 2 sYM [41]. While there have been a number of recent developments in constructing manifestly color-dual Feynman rules [42][43][44][45][46][47][48], a precise definition the kinematic algebra off-shell remains elusive. Considering the known tree-level double-copy relationships between NLSM and DBIVA and higher-derivative corrections [49,50], all of the amplitudes we compute should in principle participate in a double-copy formulation of our results. At the very least, the amplitudes we compute will serve as important checks on future breakthroughs in multi-loop studies of color-dual representations. In addition to traditional double-copy construction, we find that the quantum effective actions generated by loop effects permit a rather compact construction in terms of symmetric-structure double-copy, which was introduced in recent work by the authors [51].
3. Finally, DBIVA effective field theories touch a wide range research areas at the forefront of high energy physics. Of more formal interest, is the presence of quantum anomalies that violate the U (1) duality invariance, which the theory enjoys at tree-level [52]. The presence of these U (1) anomalies in gravity is closely linked to ultraviolet divergences computed at loop-level [53][54][55][56]. It has been argued that cancelling these anomalous matrix elements with R n counterterms can lead to enhanced UV cancellations [57][58][59]. This relationship has been demonstrated both for pure Einstein-Hilbert gravity [60][61][62], and also for less than maximal N ≤ 4 supergravity [63]. While DBIVA theories are themselves are ultraviolet divergent in D = 4, due to their simplicity they serve as an essential laboratory for probing anomaly cancellation at high loop order. Moreover, in addition to their formal theory relevance, DBI theories have garnered wide phenomenological interest in cosmology, both for sourcing nongaussianities in CMB bispectrum [64][65][66], and for their compatibility with the observed CMB tensor mode suppression [67][68][69][70][71]. Thus, with inflationary data on the horizon, understanding the perturbative structure of DBI could be particularly relevant for modeling early universe quantum fluctuations.
The outline of the paper is as follows: in section 2, we will provide a review of the onshell methods and integration techniques needed to probe the multi-loop physics studied in this paper. Then in section 3, we introduce the method of Even-point Multi-loop Unitarity (EMU) and compute the requisite two-loop tensor integrals. In section 4, we present our results, beginning with a warm-up calculation in section 4.1 where we construct integrands for NLSM through two-loop order, and compute the fully integrated amplitudes. We show that in D = 2 − 2ϵ, the leading IR divergences of the CP 1 model exponentiate. With the scaffolding of D-dimensional integration in hand, we then compute two-loop amplitudes for N = 4 DBIVA theory and pure-photon Born-Infeld theory in section 4.2. There we compute the anomalies present beyond one-loop order. Then, using the NLSM integrands of section 4.1, we perform a multi-loop double-copy to N = 4 DBIVA observables with the color-dual basis numerators available in the literature in section 4.3. After computing the catalog of twoloop amplitudes, we study the construction of quantum effective actions in section 5 as they relate to the anomalous matrix elements present in two-loop Born-Infeld amplitudes. We demonstrate in section 5.1 that the one-minus anomaly at two-loop requires the introduction of an evanescent operator at O(α ′4 ). Then in section 5.2 we take the opportunity to discuss the application of symmetric-structure double copy to cancel these anomalies, along with their relationship to higher-spin modes as studied in a recent work by the authors. Finally, in section 5.3 we lay out a Hilbert series framework for counting evanescent operators at general orders in mass-dimension. To conclude, in section 6 we discuss many directions of future work and summarize the insights gained from this study.

Color-dressed and ordered amplitudes
Here we provide an overview of the amplitudes nomenclature and organizational principles we use throughout the work. When working with scattering amplitudes, A, in an on-shell framework it is convenient to introduce the following graphical description at general multiplicity, n, and loop order, L, (2.1) As written above, S g are the internal symmetry factors for a Feynman graph, g, the propagator structure is captured by the denominator function, D g , and the theory-dependent interaction vertices determine numerator functions, N g . As written, the numerator functions carry all non-propagator kinematic information and any color-data that would be encoded in the interaction Feynman rules of the theory. By color-data we mean the typical weighting of feynman diagrams by the representation of the particles described. Most of the theories we describe are defined in the adjoint, so the color-data involves the usual dressing of vertices with antisymmetric structure constants, f abc . Although, as we will see, it will be useful to describe certain counterterms in terms of symmetric d abc color-weights, For gauge theories, it is often useful to further decompose the full, or color-dressed, amplitudes of eq. (2.1) into purely kinematic building blocks where the color information is stripped away. For example, at tree level, gauge theory n-point amplitudes can be re-expressed in terms of a trace basis decomposition as follows: A n,tree = σ∈S n−1 Tr(T 1 T σ(2) · · · T σ(n) )A n (1, σ(2), ..., σ(n)) , (2.2) where the sum is taken over the inequivalent orderings of group theory generator traces. A similar decomposition applies to one-loop gauge theory amplitudes, which we will describe in section 4. The color-ordered functions, A n , are called partial amplitudes, as they only contain on-shell information for a particular color ordering. However, like full amplitudes, they factor on poles to lower-point partial amplitudes. For vector theories partial amplitudes are independently gauge invariant since they weight linearly independent color traces. If we consider the kinematic gauge-invariant ordered amplitudes themselves, many gauge theories reveal hidden redundancy. The (n − 1)! partial amplitudes for a large catalog [26] of adjoint gauge theories, including Yang-Mills and NLSM, are linearly related to a smaller set of basis amplitudes. The first such set of partial amplitude relations identified by Kleiss and Kluijf [72], known formally as KK relations, relates amplitudes different color orderings and signature to a basis of (n − 2)! partial amplitudes: where β T denotes the inverse ordering of β, and ¡ is the shuffle product that runs over ordered permutations of α and β T . More recently, Bern, Johansson and one of the authors, (BCJ) further identified a set of kinematic relations that reduced the size of n-point partial amplitude basis down to (n − 3)! via BCJ relations [5]: Underlying both sets of amplitude relations are hidden graphical principles that manifest the amplitude level redundancy of eq. (2.3) and eq. (2.4), which we describe in the next section. Before moving on, we emphasize that everything described thus far applies in arbitrary spacetime dimensions. As such, these building-blocks are well suited for constructing loop integrands compatible with dimensional regularization.

Color-Kinematics Duality and the Double-Copy
To apprehend this graph-based structure for color-dressed theories, we first rewrite the full amplitude of eq. (2.1) as a sum over cubic graphs, as follows: 5) where N g and C g are the kinematic numerators and color factors, respectively, and Γ (3) n,L denotes the set of cubic diagrams at n-point L-loop. In this formulation of the amplitude, contact diagrams are absorbed into cubic graphs by multiplying by inverse propagators. The procedure of assigning contact diagrams to cubic graphs is referred to as generalized gauge freedom [5]. In practice, there are many different ways to assign contact diagrams to cubic graphs.
The gauge theories that we study in this paper have sufficient generalized gauge freedom to assign contact diagrams to kinematic numerators such that at every multiplicity at tree level one can write the amplitude in terms of N g that obey the same algebraic relations as the color factors, C g . For color factors that are composed of adjoint structure constants, f abc they must be anti-symmetric around vertices, and be related by Jacobi relations about each edge. For SU (N c ) we can take f abc ∝ Tr(T a [T b , T c ]) for arbitrary representation T of the gauge group.
We will also consider generalization beyond the adjoint to color-weights that can include d abc and associated algebraic relations. For SU (N c ) we can take Symmetric-structure color factors obey symmetric algebraic relations rather than anti-symmetric ones.
When a gauge theory can be expressed in a form such that the color and kinematic factors obey the same algebraic relations, we call such theories color-dual. For ordered amplitudes, the ability to find kinematic numerators that obey antisymmetry and Jacobi relations is one-to-one with the ordered amplitudes satisfying KK and BCJ relations, respectively.
An important consequence of the duality between color and kinematics is the ability to double-copy color-dual numerators with each other [5], while simultaneously preserving factorization and gauge invariance of scattering amplitudes. Because cubic-graph color-weights are not linearly independent, gauge-invariance of the full amplitude is encoded in algebraic relation between the color factors. Thus, equipped with a set of color-dual numerators,Ñ g , that obey the same algebraic relations as C g , we can make the simple replacement, C g →Ñ g , giving a new gauge invariant amplitude, M, as follows: Note that M is colorless. It satisfies manifest gauge invariance on the N g side of the double copy due to the algebraic properties ofÑ g . In the event thatÑ g also belongs to a vector theory, then M describes amplitudes in a gravitational theory. In such a case, the gauge invariance of constituentÃ and A in eq. (2.5) conspire to generate linearized diffeomorphism invariance in eq. (2.6). As long as both kinematic factors obey the same algebraic constraints as the color algebra, this double copy construction works for integrands at the multi-loop level since it globally manifests color-kinematics on all possible unitarity cuts [5].
Throughout the text we will reserve M to denote double-copy amplitudes. As shorthand, when the kinematic weights of two theories, X and Y , participate in the double copy construction of M XY we will use the outer-product to mean double copy construction as in, (2.7) At tree-level, such double-copy construction in the adjoint case can be understood equivalently [73][74][75][76][77][78] in terms of a KLT or momentum kernel matrix, S(a|b), and a BCJ spanning set of ordered amplitudes for each of the X and Y theories, The sum runs over both sets, a and b, of (n − 3)! distinct ordered amplitudes. No such momentum kernel has yet been constructed for the symmetric case [51] . As we will see, double-copy construction between symmetric color-dual kinematic numerators [51] is quite natural for capturing Born-Infeld counterterms needed for anomaly cancellation beyond one loop. Indeed, many of the higher derivative counterterms captured by double-copying symmetric kinematic factors are inaccessible 1 to the traditional double copy between adjoint kinematics. We will describe this non-adjoint color-dual double-copy in more detail in section 5.2, where we will explicate the tension between local symmetric numerators and higher-spin adjoint numerators.

4D Spinor Helicity vs. D-dimensions
All the amplitude methods we have discussed thus far are valid in arbitrary dimension. This has a number of advantages that we will touch upon when reviewing integrand construction towards dimensionally regulated amplitudes in section 2.6. While most of our calculations will be carried out completely agnostic to the spacetime dimension, the vector amplitudes we compute will have 4D symmetries that hidden in D-dimensional kinematics. It is useful therefore to also consider explicitly four-dimensional kinematics.
When working with bosonic kinematics in arbitrary dimensions, we will employ formal Lorentz covariant polarizations, ε µ a , and momenta, k µ a . The mostly minus signature will be used throughout. The only restrictions we will place on physical momenta and polarizations are the standard on-shell constraints: momentum conservation and the null-momenta condition, Furthermore, we will take external polarizations ε a formal and transverse throughout: Dimensions will generically be away from four, in D = 4 − 2ϵ, with ϵ arbitrary. As noted, it will be convenient to take D → 4 in certain circumstances after integration. In such cases we will employ appropriate spinor-helicity variables. Here we apply the same conventions of Ref. [79], which we quote now. For massless momenta, k a and k b , we have the component definition of our spinor bracket that are consistent with the conventions above, where the a i are component values of the four-vector, k µ a = (a 0 , a 1 , a 2 , a 3 ). Four-dimensional polarization dot products with fixed helicity states can be mapped as follows: where q µ is some null reference momentum that all polarizations are projected along. While all four-dimensional equivalence relayed in this paper can be achieved analytically, in practice it is often much more convenient to verify equivalence numerically.

One-loop integral basis and tensor reduction
It is well known that one-loop amplitudes are spanned by a small basis of scalar integrals [4,80,81]. This property can be exposed via a D-dimensional integral reduction algorithm due to Passarino and Veltman [82]. At one-loop, the basis of irreducible scalar products (ISPs) that include loop momenta grows in lock step with the number of propagators. Explicitly, an N -gon integral can have N − 1 factors of (k i · l), due to momentum conservation, and exactly one factor of l 2 ; this matches the number of N -gon propagators. Furthermore, factors of (ε i · l) can be mapped to tensor structures of external kinematics with factors of (ε i ·k j )(k i ·l)/(k i ·k j ) using polarization completeness relations [57] . Thus, any N -gon tensor integral permits a partial fraction decomposition in terms of inverse propagators, and thus can be mapped to a scalar basis of integrals: where T µ 1 µ 2 ...µr and C M are functions exclusively of external kinematics. This relationship holds D-dimensionally for any r-rank N -gon one-loop integral. Thus, any one-loop amplitude can be expressed completely in terms of a basis of one-loop scalar integrals. This special property of one-loop integration is not the case for generic loop order, where the basis of ISPs grows faster than number of inverse propagators. As such, any universal integral basis beyond one-loop must include integrals with non-unit powers of the propagators. The remarkable simplicity of one-loop amplitudes is even more dramatic in a fixed spacetime dimension. In D = 4, the kinematic projection to a basis of scalar integrals saturates at the box integral due to the appearance of Gramm determinants. Roughly speaking, N -gon integral coefficients of eq. (2.15) come dressed with a factor of G N , where k i are the external momenta flowing into the N -gon integral. In fixed spacetime dimension, D, the momenta are D-component vectors and thus the Lorentz product matrix (k i · k j ) must have a null space for N > D. As a result, C M ̸ = 0, only when M ≤ D + 1. Furthermore, in D = 4 − 2ϵ dimensions, any scalar pentagon integral can be rewritten up to O(ϵ) in terms of five pinched scalar box integrals [83]. All the information of one-loop 4D amplitudes can therefore be completely determined by evaluating up to box integrals, regardless of the multiplicity. We will demonstrate in the next section how this can be used in 4D integral construction.
The last important feature of one-loop basis integrals is that many of them can be evaluated in arbitrary dimension for a very general set of parameters appearing in the exponents of the propagators. Take for example the D-dimensional massless triangle and bubble integrals [84] that we will use throughout the text: where we have introduced the subscript notation for indexed sums X 12...n = X 1 +X 2 +· · ·+X n . Throughout the text we will suppress the regularization scale that appears in the argument of the logarithms as they can be restored by dimensional analysis. The above integral expressions hold for massless external kinematics, K 2 i = 0, which applies to all the integrals needed for the physical processes described in the text.
The existence of closed form expressions, like those in eq. (2.17), while common for oneloop integrals, are incredibly rare for higher loop order outside of very specialized kinematic regimes. Luckily, the integral basis at four-point two-loop order for even-point theories can be reconstructed from the one-loop expressions in eq. (2.17). This is a special property of what we refer to as recursively one-loop amplitudes. We can exemplify this property in the simple context of scalar φ 2k theory.
Consider the interactions needed to construct two loop amplitudes for the generic arbitrarily weighted φ 2k -scalar theory: The two-loop amplitude for φ 2k theory is simply generated by evaluating the following set of scalar 1-particle-reducible (1PI) integrals: Since the external states are massless, only the first two integrals, those appearing in Fig. 1, are non-vanishing in dimensional regularization -the rest are scaleless. Therefore the above even-point perturbative predictions of φ 2k are equivalent to φ 4 theory at four-point through two-loop order: (2.20) In addition to being perturbatively equivalent to the quartic sector of the theory at two-loop, we can see above that the integrals for the Feynman diagrams in eq. (2.19) are recursively one-loop. In other words, to integrate the full amplitude at two-loop, we only need the iterated basis of one-loop integrals. This observation comes with the added advantage of permitting the application of one-loop Passarino-Veltman tensor reduction [82] to two loop integrals. This property of recursively one-loop integrals will dramatically simplify the higher derivative even-point theories that we consider.

Even-point Effective Field Theories
In this section we provide a catalog of the even-point effective field theories studied throughout the paper. Some of the theories we consider, like Volkov-Akulov and supersymmetric extensions of Born-Infeld, have natural formulations in terms of veirbeins that depend on fermionic degrees of freedom. For such theories, it is useful to distinguish between the coordinate independent actions, S, and the coordinate-dependent Lagrangian densities, L, Nonlinear Sigma Model The simplest effective field theory that we will consider is the nonlinear sigma model (NLSM). The chiral NLSM Lagrangian density with diagonal symmetry group G = SU (N c ) can be expressed in terms of the chiral current, j µ = U † ∂ µ U , as follows, where the trace is over the color indices of the gauge group, π ≡ π a T a . On the right hand side we have applied the Cayley parameterization to the SU (N c ) group elements, U = 1+π/fπ 1−π/fπ , where f π is the dimensionful pion decay width. By construction, this model has a left-right global symmetry, G×G, that is nonlinearly realized, and a linearly realized diagonal subgroup, G. This theory is the unique two-derivative EFT that is invariant under constant shifts of the pion field in color space: π a → π a + c a . (2.23) We could also consider higher derivative generalizations of the NLSM that are invariant under the shift symmetry of eq. (2.23). Such EFTs correspond to the so-called chiral limit of chiral perturbation theory [85][86][87][88], describing massless pion scattering below the chiral symmetry breaking scale of QCD: (2.24) By construction, each of the higher derivative operators appearing a bove are invariant under the shift symmetry of eq. (2.23). These higher derivative generalizations of NLSM have been studied both in the context of the S-matrix bootstrap [89][90][91], and soft theorem bootstraps of Refs. [8][9][10][11][12]. Given the relevance of eq. (2.24) in Standard Model pion scattering, the amplitudes for G = SU (2) isospin pions have been computed through two-loop order [92][93][94]. In this work, we add to the literature by computing the two-loop amplitudes of eq. (2.22) for arbitrary SU (N c ) color structure. We then use these amplitudes to compute N = 4 DBIVA amplitudes via the double copy in section 4.3.
In addition to the chiral NLSM described above, our methods apply equally well to a target space model describing pion dynamics on a coset manifold G/H, with a nonlinearly realized global symmetry G broken down to a linearly realized isometry group H. Recent literature that study the universal soft behavior of these models can be found in Refs. [95][96][97][98]. One such model that we will study in the text is the CP N model where the global SU (N + 1) is spontaneously broken down to a local U (N ) symmetry on the target space. The Lagrangian for the theory we provide below: where z = (z 1 , z 2 , ..., z N ) is a complex vector in the fundamental representation of U (N ), and P (z,z) is the Fubini-Study metric on CP N . The full SU (N + 1) symmetry of the theory is realized due to the imbedding of CP N in a complex space one dimension higher, C N +1 . In section 4.1, we will demonstrate through explicit calculation that the leading IR divergences of the CP 1 model exponentiates in D = 2, the critical dimension of the theory.

Born-Infeld Photons and Dirac Scalars
The second even-point EFT we will study is Born-Infeld theory [99]. This theory describes the dynamics of open string endpoints with U (1) charge attached to a D-dimensional spacetime with Dirichlet boundary conditions. Given the stringy dynamics orthogonal to the spacetime, the electric field generated by the charged endpoints has a maximum field strength of order the inverse string tension, 1/α ′ . The effective Lagrangian for this theory A close cousin of Born-Infeld theory is the dimensionally reduction to Dirac scalars, which describes a D-dimensional spacetime propagating in a (D + 1)-dimensional background: where φ is the spacetime coordinate in the orthogonal space. This theory is uniquely determined by considering the most general polynomial of X = (∂φ) 2 , P (X), constrained to be invariant under the field redefinition [10]: When the transformation law leaves the polynomial P (X) theory invariant, there is an additional set of D + 1 conserved charges, and the Poincaré symmetry of P (X) EFT is promoted from D to (D + 1) dimensions. One can verify that applying eq. (2.28) to eq. (2.27) shifts the Lagrangian by a total derivative [100]. Similar to the soft theorem constraints on the building blocks of NLSM and χPT, the DBI action can be bootstrapped by requiring the amplitudes vanish at O(p 2 ) when taking external momentum, p, soft [10].
Volkov-Akulov Fermions and Supersymmetry Now we will describe the supersymmetric extensions of the Born-Infeld action given above in eq. (2.26). First we start by defining the theory of shift symmetric fermions first written down by Volkov, Akulov (VA) [101].
The easiest way to define this shift symmetric VA theory is by constructing a volume form in the veirbein frame out of manifestly supersymmetric 1-forms. That is, define a 1-form, ω m = e m µ dx µ , where the veirbeins (frame fields), e m µ , are expressed in terms of the fermion fields, λ, as follows: This frame field is manifestly invariant under the superspace shift: The Volkov-Akulov action in the veirbein frame is then simply the fermionic world-volume integral over the D-form [102], This is not dissimilar from the construction of the pure scalar DBI action above. Likewise, we can write the Lagrangian in the spacetime coordinate frame in parallel to eq. (2.26) and eq. (2.27). Transforming back into spacetime coordinates yields the following Lagrangian density for self interaction VA fermions upto four-point interactions: Furthermore, a convenient choice of field redefinition was proposed by Komargodski and Seiberg (KS) [103,104], which defines a new fermion field, ψ = ψ(λ,λ), in a way that is perturbatively equivalent to the VA fermions. This construction is equivalent to a theory of nilpotent chiral superfields [105][106][107], relevant for the construction of α-attractor models of inflation [108]. Applying the KS nonlinear field redefinition yields the considerably simpler Lagrangian for VA theory: With the world-volume formulation of VA theory in hand, we can also trivially construct the action for maximal N = 1 supersymmetric Born-Infeld theory in D = 10 [109]. We simply replace the flat space background of eq. (2.26) with the fermionic background of eq. (2.32), = ω 1 ∧ ω 2 ∧ · · · ∧ ω 10 det(η mn + αF mn ) .

(2.34)
It is worth noting that the fermion-vector interaction is introduced via the non-minimal coupling between the field strength and the fermionic spacetime metric [110]. This differs from the linearly realized supersymmetry of super-Yang-Mills theory, which minimally couples the vector mode to the chiral fermion. This non-minimal coupling is a requirement of supersymmtry due to the even-point nature of abelian vector theories like Born-Infeld. The action above in eq. (2.34) can be similarly re-expressed in terms of our canonical spacetime coordinates, giving us the following flat-space Lagrangian for maximally supersymmetric Born-Infeld theory in D = 10: (2.35) To recover the D = 4 theory, we simply dimensionally reduce in a way that preserves the maximum number of supercharges, which is N = 4 in D = 4. Implementing a dimensional reduction that preserves the supercharges of eq. (2.35) is rather straightforward in an on-shell framework: Since the double copy construction holds D-dimensionally at tree-level [26], double-copying maximal sYM with NLSM will yield the analogous spectrum of maximally supersymmetry DBIVA theory, whether in D = 10, as described above, and in D = 4.
Abelian Open String Similar to the amplitudes of DBIVA theories described above, the tree-level amplitudes for the open superstring (OSS) [111,112] likewise permit a double copy construction [113]. To construct open-superstring amplitudes we simply double copy maximally supersymmetric super-Yang-Mills (sYM) and Chan-Paton dressed Z-theory amplitudes [50,114,115], as follows, where ⊗ implements the field theory double copy [5,116] defined in eq. (2.7). Moreover, the amplitudes of DBIVA emerge from the above construction as the field theory limit of the abelian open superstring [117]. To recover the abelian sector, one simply sums over all Chan-Paton color orderings on the side of the bi-colored Z-theory amplitudes that obey string monodromy relations. In the field theory limit, the observables for so-called abelian Z-theory are simply the amplitudes generated by the NLSM Lagrangian [50]. That is, given a n-point amplitude in Z-theory, with field theory ordering a = (a 1 , a 2 , ..., a n ) and string theoretic color ordering A = (A 1 , A 2 , ..., A n ), one finds the following relation: A NLSM (a 1 , a 2 , ..., a n ) ≡ lim where Z × are abelian Z-theory amplitudes of [50], defined by summing all possible orderings of Chan-Paton factors: ..,A n−1 ,n) (a 1 , a 2 , ..., a n ) . (2.39) Since the field theory limit of abelian Z-theory produces NLSM amplitudes, the field theory limit of the abelian OSS, gives rise to DBIVA observables at leading order in α ′ : Additional details on the bi-colored Z-theory amplitudes these string-theoretic double-copy constructions can be found in [50,[113][114][115]118] and references therein. We emphasize that while NLSM and DBIVA are the two theories that we will focus on, the methods that we develop in the next section apply much more broadly to the full catalog of even-point effective field theories.

Duality Invariance and Higher Derivatives
The four-dimensional amplitudes of DBIVA theories have an additional special property that we will study in his work. At tree-level, the amplitudes generated by eq. (2.26) and eq. (2.34) exhibit U (1) duality invariance [119,120]. A theory is said to exhibit duality invariance when the effect of rotating field strengths, F µν , into the dual-fields, G µν , defined below, leaves the equations of motion invariant [121,122]. This dynamical symmetry gives rise to a 4D helicity selection rule [120] whereby tree-level matrix elements vanish on-shell outside the aligned-helicity sector, where n (+) and n (−) is the number of external positive and negative helicity photons, respectively. For N = 0 pure Born-Infeld theory [99,123], this symmetry is violated at four-point one-loop in the all-plus (+ + ++) helicity sector [124]. In the text, we will show how this anomaly is propagated to two-loop, and demonstrate that U (1) duality invariance is further broken in the (− + ++) configuration beyond one-loop. Before proceeding, we comment briefly on the higher derivative corrections one might expect to appear in duality invariant effective field theories. In the case of supersymmetric DBIVA at four-points, U (1) symmetry is promoted to an R-symmetry, and is thus protected perturbatively to all orders by supersymmetric Ward Identities [125]. We show this explicitly in the case of N = 4 DBIVA through two-loop in section 4.2. However, this supersymmetric enhancement of U (1) duality does not in general apply to higher multiplicity amplitudes.
While the leading order in α ′ must be duality invariant at all multiplicity for DBIVA [21], beyond the leading order, R-symmetry permits non-vanishing matrix elements outside the split helicity sector. The low energy effective action of the OSS is one such exemplar of an EFT that respects U (1) duality at leading order in α ′ , but produces duality violating matrix elements at higher-orders above four-point [50].

On-shell Unitarity methods
Above we have reviewed the known behavior of one-loop integral bases in D = 4 spacetime dimensions. In this section we will briefly demonstrate why traditional 4D unitarity methods are so powerful -and why we must go beyond them to accommodate the multi-loop calculations of interest here.
Integrands from 4D Cuts Applying one-loop integral reduction in 4D tells us that any one-loop massless amplitude, for which tadpoles are scaleless, can be expressed in terms of box, triangle, and bubble diagrams, with each integral weighted by dimensionally-dependent functions of kinematics: where the scalar basis-integrals, I D N , are evaluated in D = 4 − 2ϵ, and R is a rational term that emerges from series expanding around ϵ = (4 − D)/2. To construct integral coefficients, C N , we make use of the Optical Theorem for the S-matrix, which states that, where S = 1 + iT . Applying this relation to eq. (2.43) allows us to identify the imaginary part of A 1-loop with the imaginary parts of each basis integral weighted by the respective 4D coefficient, C N (4): Using the optical theorem to constrain scattering amplitudes is known as the unitarity method [1,3]. In order to isolate each C N (4) appearing in eq. (2.43), we must place particular internal legs on-shell. That is, for each internal loop propagator that we wish to take on-shell, 1/P 2 i , we make the following replacement inside the integral: By replacing internal propagators with positive-root delta functions, one can isolate the imaginary parts needed to determine the integral coefficients. This approach of taking individual legs on-shell is known as the method of generalized unitarity -and it can be systematically applied to construct all one-loop amplitudes directly from 4D unitarity cuts [4]. Using generalized unitarity, the integral coefficients correspond to the following set of iterated cuts: , (2.49) where exposed legs are summed over all internal states crossing the cut, and the labels P i are external momentum scales flowing into the specified cuts. The subtractions account for the overlapping information in the respective N -gon cuts. By momentum conservation, the momentum flow in the unlabeled vertex is completely determined by the other N − 1 vertices.
As a concrete example, consider the 4D cut construction of the bubble coefficient for pure Born-Infeld theory: where in the last step we have made the replacement (l 1 + l 2 ) = −(k 1 + k 2 ) by momentum conservation. Determining the leading behavior of the 1-loop Born-Infled amplitude thus amounts to specifying all possible 4D cuts. While this method is incredible powerful, as we can see from eq. (2.45) it misses a key piece of the amplitude -namely, the rational terms, R. For some theories like N = 4 super-Yang-Mills in D = 4 [2], rational terms are absent because the loop power counting for an N -gon cut is bounded to be less than N − 2. However, for generic theories of interest in this paper, like Born-Infeld, rational terms are present and necessary for extracting anomalous matrix elements and higher-loop sub-divergences. While there are methods developed to extract these rational terms at one-loop using µ-integrals [80], at generic loop order evaluating these µ-integrals can become fairly complicated. For this reason, to perform the multi-loop calculations in this paper, we will employ cut construction in general dimensions [126,127].
Integrands from D-dimensional Cuts The presence of 4D rational terms in the amplitude is due to the appearance of (D −4) factors. By constructing integrands D-dimensionally, we can guarantee that our integral coefficients are sensitive to these overall factors. Rather than summing over internal 4D states, we will instead compute D-dimensional cuts using the formal polarizations described in section 2.3. The D-dimensional analog of eq. (2.50) asserts that the one-loop integrand I 1-loop of pure Born-Infeld theory should satisfy the following constraint: where Cut(I 1-loop ) extracts the kinematic numerator by imposing a maximal cut on the oneloop integrand . (2.55) To carry out the sum over formal polarization states in eq. (2.54) we apply the following D-dimensional completeness relation for on-shell polarizations, with null reference momentum, q 2 = 0. The dependence of q is a gauge choice that disappears when on-shell kinematic constraints are imposed. At loop level, the state-sum of eq. (2.56) can potentially lead to dimension dependence in the D-dimensional integrand. This is due to terms that contain products of internal polarizations: where D s is the spin-dimension of the internal photons, and thus D s −2 just counts the number of photon states. To calculate our loop-level amplitudes, we use dimensional regularization to regulate the integrals. We will adopt the conventions of [128,129] and take the spin-dimension, D s , to be the same as the dimension appearing in dimensional regularization, D = 4 − 2ϵ.
We use them interchangeably throughout. The appearance of these D-dependent factors in the state sum are the source of rational terms in the Born-Infeld S-matrix. Throughout the text, we will use the generalized unitarity cut condition of eq. (2.54) to fix our multiloop integrands. With this in hand, in the spirit of Forde's one-loop cut construction, we will now systematize the procedure for capturing all the cut information needed for even-point theories at the multiloop order.

Even-point Multi-loop Unitarity
In this section we will describe the methods we have developed to compute the multi-loop results of this paper. We begin by introducing the notion of Even-point Multi-loop Unitarity (EMU), which is an organizational principle at the foundation of our unitarity-based integrand construction. EMU is an extension of the method of maximal cuts [130], which is a hierarchical approach to perturbative calculations [131]. EMU is a constructive algorithm aimed at capturing all the perturbative information needed at general loop order and multiplicity. We describe the algorithm below, and provide a 2-loop 4-point example that we choose to study in this paper.

Even-point Multi-loop Unitarity (EMU)
• Step 0 -At n-point l-loop, enumerate all 1-particle-irreducible (1PI) graphs constructed from (n + 2l − 2)/2 four-point vertices. We call these diagrams the maximal cut (MC) diagrams. Graphs with higher-point blobs are grouped into the N k MC category, where k is the number of collapsed internal propagators relative to the MC graphs.
• Step 1 -Culling from N k MC → ∂N k MC: -Step 1A -Discard all diagrams at the given order N k MC that capture scaleless behavior from internal kinematics. For massless theories, this amounts to throwing out diagrams that contain one of the following internal nodes: Step 1B -Discard all diagrams at the given order N k MC that capture scaleless behavior from external kinematics. At n-point, this amounts to rejecting diagrams that contain n − 1 external edges attached to a single blob: If the culling procedure of Step 1 produces an empty set of graphs, ∂N k MC = ∅, then the routine terminates. If not, collapse one of the internal propagators for all diagrams in all topologically distinct ways. Collapsing a propagator simply means merging the two nodes connected by that propagator's edge and removing that edge all together. So the resulting graph will have one less internal edge and one less internal node. This collapsing step will take the ∂N k MC set of diagrams to N k+1 MC. Once there, repeat Step 1A and 1B until the routine terminates.
This procedure for collecting all the cut information using EMU can be represented diagrammatically as a sequence of culling and collapsing as follows: After applying EMU and collecting all the diagrams up to some order N k MC (l,n) , the set of cuts we use to fully constrain the n-point l-loop integrand, Ω (l,n) , is the intersection of these sequence of unitarity cuts: Taking the intersection accounts for the overlapping information stored in the cut diagrams 2 at different orders in N k MC. As a concrete example, let's consider the case of two-loop fourpoint of interest in this paper. First, we obtain the following set of diagrams at the MC level built from three four-point vertices: Each blob corresponds to an on-shell tree-level amplitude from the even-point theory of interest. The third diagram is scaleless on support of dimensional regularization, and thus we discard it in Step 1A. This gives the following restricted set of MC diagrams: Since ∂MC (2,4) is not empty, we proceed to Step 2. After throwing out the scaleless diagram in MC (2,4) , this leads to the following set of N 1 MC diagrams for two-loop four-point: Now the first diagram is discarded in Step 1A, and the second diagram is discarded in Step 1B. Thus, we need not consider any next-to-maximal cuts for two-loop four-point amplitudes in even-point theories. Since ∂N 1 MC (2,4) = ∅, this concludes the EMU cut construction.
As we can see from the diagrams that contribute at two-loops four-point in eq. (3.6), evenpoint theories can be fully constructed from convolutions of one-loop integrals. As described earlier, we call such integrals recursively one-loop since we can iteratively apply one-loop unitarity and tensor reduction methods in order to recover the full multi-loop structure. In the next two sections, we first illustrate the recursive behavior of convolution integrals, and then use their properties to develop multi-loop tensor reduction methods that we use throughout the paper.

Multi-loop recursive integrals
As we found above in our example of EMU, the only diagrams that contribute at two-loop four-point are those given in Fig. 1. Due to the bubble integral insertions, both of these integrals can be evaluated in terms of the one-loop basis integrals of eq. (2.17). This is a direct consequence of the single-scale nature of bubble insertions. Since internal bubbles are single-scale, they will only contribute additional powers of inverse propagators. This has the effect of shifting an L-loop convolution integral to (L − 1)-loop order, with a propagator power determined by the mass-dimension of the evaluated bubble integral. Explicitly, we can evaluate a D-dimensional nested bubble with internal momentum, q µ , flowing into the diagram, where α 1 and α 2 are the degree of the denominators appearing in the bubble integral. For example, consider a banana integral diagram constructed from six-point contacts that would appear in a three-loop application of EMU with external momentum scale q µ = (k 1 + k 2 ) µ .
Evaluating the simplest case where all α i = 1, one finds that recursively applying eq. (3.8) yields the following sequence: For simplicity, we have dropped the Γ-function dependent factors at each step in the evaluation. This procedure of evaluating bubble integrals at successive orders also works for tensor bubble-integral insertions as well. As we'll show in this next section, such tensor bubble insertions can be reduced to scalars integrals via Passarino-Veltman [82]. While there are a number of tensor reduction algorithms available in the literature [132][133][134][135][136][137][138][139], our ability to use Passarino-Veltman on the recursive integrals of Fig. 1 will prove useful for evaluating terms with factors of (ε i · l j ), which include formal polarizations. This procedure will aid in our pursuit of extracting the rational part of the two-loop Born-Infeld amplitudes that we study in section 4.2.2.

Two-loop tensor reduction
First we will study the canonical one-loop case needed for the nested bubble integrals described above, and then show that it naturally generalizes to the two-loop for the recursive integrals of interest. We begin with a general expression for the rank-n tensor bubble: where we have introduced the following shorthand for the symmeterized rank-(m+2k) tensor: By contracting the tensor integral with the internal scale, K µ , and the metric, η µν , we obtain the following two constraint equations: Performing the same contractions on the symmeterized tensor defined above in eq. (3.12) yields the following: (3.16) Figure 2. Graphical depiction of the the double-bubble integral. Every exposed internal edge represents a propagator in the integrand.
These linear relations can be used to line up the coefficients with distinct tensor structures, which gives the following: 18) and similarly so for the metric contraction: Treating the distinct tensor structures as basis elements, we thus conclude the following set of linear relations between the coefficients for the rank-n tensor bubble: where D = η µν η µν . These constraints can be rearranged to give the following recursive definition for the coefficients a (m,k) : The base step is simply the scalar bubble integral, a (0,0) = I 2 (K). Constructing a (m,k) from this recursive definition, the one-loop integrand of eq. (2.54) that contains factors of (ε i ·l) and (k i ·l) can be expressed completely in terms of the bubble integral, I 2 (K), weighted by dimension dependent numerical factors and external vector kinematics. We note that these coefficients are also sufficient to evaluate the contribution from the double-bubble integral whose propagator structure is sketched in Fig. 2. This is best demonstrated with an explicit example. Consider the following integral, I ex.
2bub , that functionally captures terms that could appear in the cut construction of the double-bubble integral at two-loop: where v i is a stand-in for external kinematics, ε i or k i . While the numerator mixes factors of l 1 and l 2 , the denominator can be separated. This allows us to re-express the above integral completely in terms of iterated tensor bubbles of eq. (2.17): Then, by applying eq. (3.11), and plugging in the expressions for a (m,k) , the kinematic numerator of I ex. 2bub no longer mixes loop momenta. Thus, the integral is separable and can be expressed as a product of scalar bubbles integrated over l 1 and l 2 .
A similar procedure can be applied to the ostrich type diagrams of Fig. 3. However, the integration procedure is a little more delicate than the I ex.
2bub could be expressed as an iterated bubble, we only needed to consider integer powers of the denominators when constructing the recursion relations. In contrast, ostrich integrals will lead to non-integer, ϵ-dependent powers of loop propagators. We are therefore interested in performing an x-dependent tensor reduction on the triangle integral, with K 12 = K 1 + K 2 introduced as shorthand notation. We note that x is a non-integer value that will inherit dependence on ϵ from integrating over the internal bubble, The symmetrized triangle tensor takes the following definition: Given this definition, when we perform the tensor contractions over K 1 and K 2 , the degree of x will get shifted: For our purposes, external momenta are taken to be null, K 2 1 = K 2 2 = 0. By construction, contracting with the metric will yield a scaleless integral, giving us the following constraint: Note that in the above contractions with K 1 and K 2 the degree of the denominator get's shifted from x → x − 1. For normal one-loop triangle integral reductions for which x = 1, this would lead to a tensor bubble that we have computed in the previous section, However, since the two-loop integration leads to ϵ-dependence, this will now induce additional factors of scalar I 3,ϵ−n final tensor reduction. Proceeding with the computation, these integral constraints can be similarly applied to our symmeterized tensors, yielding the following set of contractions: Keeping this in mind, we obtain the following linear relations between the coefficients, (3.37) Just as was done for the scalar bubble, these functional expressions can be rearranged to construct family of recursion relations for the tensor coefficients: This system of equations uniquely fixes the rank-n massless triangle tensor integral in eq. (3.26) with non-integer powers in the denominator. Using these integral reduction formulae, both the one-loop and two-loop integrals can be expressed completely in terms of scalar one-loop bubble and scalar triangle integrals, which we are provided in eq. (2.17) of the previous section. Before proceeding to our results, we will introduce a bit of notation that will be relevant for capturing loop-level contributions to the vector amplitudes.

Gauge-invariant on-shell basis
The final bit of machinery we introduce for the multi-loop calculations of this paper is a spanning set of D-dimensional four-photon on-shell tensor structures. All of the vector integrands constructed in the following sections will be fixed on D-dimensional cuts, and projected to a basis of D-dimensional photon tensors.
Being agnostic to dimension allows for a particularly algorithmic approach to identifying rational terms. In addition, by projecting to a D-dimensional basis of gauge-invariant tensor structures, we can more easily track the one-loop divergences that propagate to two-loop order. This is due to the existence of evanescent operators, those which vanish when plugging in 4D states, but are non-vanishing in general dimensions. As we will see in pure Born-Infeld amplitudes of section 4.2, these so-called evanescent operators are the cause of two-loop divergences that are obscured when looking at the one-loop behavior in D = 4 exclusively. Tracking divergences introduced by evanescent operators has been an active area of research Standard Model effective field theory (SMEFT) [140][141][142][143][144], the UV behavior of quantum gravity [57,145,146], and more generally as an area of formal theory interest [147][148][149][150][151]. Below we will give a concrete example of an evanescent operator expressed with notation used in the text to provide some justification for our D-dimensional tensor basis.
First, we will use the pair of D-dimensional tensor structure, f ij f kl and f ijkl defined previously in [51], where F µν i = k µ i ε ν i − k ν i ε µ i are linearized field strengths, and the trace is taken over spacetime indices. With these vector building blocks, there are exactly four D-dimensional four-photon matrix elements one can write down at O(k 8 ): where the Mandelstam invariants are defined as s ij = (k i + k j ) 2 . The index subscripts correspond to powers of Mandelstam invariants that are easly generalized to span arbitrarily high mass-dimension, Immediately we can see that the full D-dimensional basis must have a non-trivial null space when projected down to D = 4. Given that the 4D helicity space is overdetermined, we are able to define the following evanescent amplitude T ev. that vanishes when constrained to any 4D states: Using the helicity projection rules of eq. (2.14), one can verify that T ev. indeed vanishes in D = 4, while clearly non-vanishing in general dimension. This behavior will be particularly relevant for interpreting the loop level results where D = 4 − 2ϵ. Indeed, it is critically important that we track D-dimensional contributions when probing for divergences at multi-loop order. Consider a two-loop integral for which we need to integrate the following quantity: (1,2,l 1 ,l 2 ) T arb.
(l 1 ,l 2 ,3,4) wherel 1 = −l 1 = l andl 2 = −l 2 = −(l + k 12 ). We have dressed the darker vertex with the evanescent contribution, T ev. , which could be the result of an unspecified loop integral, and take a generic arbitrary vector structure T arb. to dress the lighter vertex. Exposed legs are taken to be on-shell in the numerator of this expression. To see what can go wrong with relying only on four-dimensional cut information, we can consider a particularly simple case. Take the vector insertion on the right hand side to be T arb.
(1234) ≡ f 12 f 34 . By applying the state sum and tensor reduction formulae from the previous sections, this can be evaluated explicitly. The result is given, where I 2 (k 12 ) is the scalar bubble integral in the s 12 -channel, and in the second line we have plugged in D s = 4 − 2ϵ. Thus far, this is exactly what we should expect. Since the amplitude T ev. vanishes in D = 4, the 4-dimensional cut above must vanish. The Optical Theorem thus disallows imaginary parts from logarithms that would appear post integration. This physical constraint is imposed by the factor of (D − 4), which absorbs the divergence of I 2 (k 12 ) and pushes the logarithm to be subleading at O(ϵ). However, suppose that T ev. came dressed with a 1/ϵ-divergence from some nested integral in a full two-loop calculation. That is, suppose the operator insertion came from the leading divergence of a one-loop integral, such that, Now the two-loop contribution would be divergent despite the vanishing 4D cut between the evanescent operator and the example operator, T f 2 , chosen above: We would have missed this if we relied solely on four-dimensional cut information.
It is straightforward to see that the three D-dimensional structures given in eq. (3.41) span the predictions of four-photon effective operators at arbitrary mass-dimension. The first two T F 2 F 2 (x,y) , and T F 4 (x,y) , span the external (± ± ±±) and (± ± ∓∓) helicity sectors. The third T F 3 (x,y) captures the (± ± ±∓) helicity configurations.
As an aside it is worth noting that the single F 3 -insertion vector permutation invariant, stA F 3 (s,t) , can be expressed concisely in terms of T F 2 F 2 (x,y) as follows, While not obvious when expressed in this form, the numerator of the above tensor structure is proportional to the permutation invariant in the denominator, T F 2 F 2 (0,2) − g 1 g 2 g 3 g 4 ∝ σ 3 , and thus stA F 3 (s,t) is local by construction. As eq. (3.41) represents a set of tensor structures which completely spans the space ofpredictions of D-dimensional permutation invariant photon effective operators, we can write an on-shell effective amplitude parameterized by numeric Wilson coefficients, a (x,y) , as follows, In Ref. [51] it was demonstrated that this photonic effective field theory (EFT) amplitude permits a double copy construction by either introducing a set of d abc -type symmetric algebraic structures, or more traditional adjoint f abc -type kinematics but at the cost of factorizing over higher spin modes. While we will comment on the double-copy properties of this amplitude in section 5.2, for now we will just stress that these D-dimensional operators will serve as an useful basis for capturing divergences and anomalies in the Born-Infeld S-matrix. With that, we are prepared to proceed to the loop level results.

Loop-level results
We now apply the EMU approach introduced in section 3 to construct two-loop amplitudes in NLSM and DBIVA theories and reduce them to an appropriate basis of integrals. We then evaluate the D-dimensional scalar integrals in the dimension of interest. We will focus primarily on D = 4−2ϵ for DBIVA amplitudes, but will also consider D = 2−2ϵ, which is the dimension for which NLSM is critical. We will then project D-dimensional tensor structures along defined 4D spinor helicity states using the conventions described in section 2.3. This will clarify where anomalous matrix elements contribute for non-supersymmetric BI, which has a classically conserved U (1) symmetry.
We begin with an instructive example of NLSM pions through the two-loop order. These amplitudes can serve as the scaffolding for constructing N = 4 DBIVA integrands using the double copy even without explicitly finding a color-dual representation. An exciting realization of our results is that the loop-level double copy construction of N = 4 DBIVA is consistent with the results that we compute directly via unitarity. This is the first application of double-copy construction for multi-loop amplitudes for non-gravitational theories. Furthermore, it suggests that there should exist a color-dual representation of two-loop NLSM integrands, a representation that has yet to be constructed explicitly.

NLSM via EMU
We begin with the color-dressed NLSM tree amplitudes that will serve as the kinematic building blocks to appear in our unitarity construction. For our purposes, a convenient color basis is the so-called half-ladder 3 basis of Del Duca, Dixon and Maltoni (DDM) [152], The half-ladder color factors are defined in terms of the antisymmetric structure constants as, In the above expression, each half-ladder color factor, associated with a particular channelgraph, dresses a color-ordered amplitude A NLSM (1σn) . At four-point, it is well known that the (ijkl) ordered amplitude for NLSM is simply, The color-dressed amplitude is of course entirely Bose-symmetric, but we use the subscript (1|23|4) to emphasize a functional choice of the (1|σ|4) half-ladder basis.

One-loop
At one-loop, the unitarity construction of gauge theory amplitudes is particularly simple, since the color factors can be decomposed into a DDM-like ordered color basis. That is, all color factors, C g associated with the following cubic representation of the integrand: can be expressed uniquely as a sum over a canonical ordering of color factors by iteratively applying the color jacobi identity. In doing so, all 1-loop color factors, C 1-loop g , can be expressed as a sum over n-gon integrand factors, C n-gon (a 1 a 2 ...an) in the following way, where β (σ) g are factors of {−1, 0, 1} depending on the color structure of the diagram, C 1-loop g . The n-gon basis diagrams take the following definition in terms of adjoint structure constants, f abc : This allows us to express the one-loop amplitudes in terms of a sum over the (n − 1)! distinct color factors, weighted by the contributions from the color-ordered Feynman diagrams: As this is a minimal basis we will consider cuts that allow a targeted identification of each color-weight's coefficients.
Construction Using the four-point tree-level NLSM amplitude with the reasoning of section 3, it is straightforward to write down an off-shell integrand associated with the one-loop bubble: We define the internal loop momenta as,l 1 = −l 1 = l andl 2 = −l 2 = −(l + k 12 ). The integrand numerator can be easily evaluated, yielding the following expression: states A NLSM (l 1 |12|l 2 ) A NLSM (l 2 |34|l 1 ) = 4f −4 π C box (1234) (k 2 ·l 1 )(k 3 · l 1 ) + (k 1 ·l 1 )(k 4 · l 1 ) (4.10) + 4f −4 π C box (1243) (k 1 ·l 1 )(k 3 · l 1 ) + (k 2 ·l 1 )(k 4 · l 1 ) . (4.11) Again, following the reasoning in section 3 we do not need to worry about any contact terms that may have been missed in this cut -any such terms do not exist for an even-point theory. We note that since the internal pions are identical, this integral comes dressed with an internal symmetry factor of S (12|34) = 1 2 . Moreover, to capture all channels we can sum over all cyclic permutations:  Now, given the appearance of loop momenta in the integrand, we must proceed to the next step in our integration procedure: tensor reduction.
Reduction Since there are up to two powers of loop momenta appearing in the integrand, it is clear that there will be dimensional dependence when applying the Passarino-Veltman integral reduction of eq. (3.23). In doing so, we can write the one-loop pion amplitude completely in terms of scalar kinematics and scalar bubble integrals: (2, 3, 4) . Integration In the study of NLSM loop-level amplitudes, we provide two examples for the integration dimension, D = 4 − 2ϵ and D = 2 − 2ϵ, since the critical dimension of NLSM is D = 2. In an MS renormalization scheme, the ϵ-expanded bubble integrals for these two dimension choices are as follows: 14) Plugging these into our D-dimensional expressions, and dropping scheme dependent rational terms at subleading order, yields the following color-ordered 1-loop pion amplitudes for each respective dimension: where as defined in section 3.3, σ 2 = (s 2 12 + s 2 13 + s 2 14 )/8, and the full amplitude is recovered by summing over cyclic permutations of (2,3,4), It is worth emphasizing that while the divergence in D = 4 − 2ϵ is an ultraviolet divergence, the one in D = 2 − 2ϵ is a logarithmic IR divergence, akin to that of N = 4 super-Yang-Mills (sYM) in the critical dimension of D = 4 − 2ϵ. Since all the particle states above are scalars, there are no helicity structures to map into, and thus we will bypass the final integration step of projection. Equipped with this one-loop warmup, we are now prepared to take on a two-loop example.

Two-loop
The two-loop calculation is exactly the same procedure as one-loop, with some additional details that need to be accounted for when performing the integration step. We begin just as above with integrand construction via D-dimensional unitarity.
Construction The color-decomposition becomes slightly more complicated at two-loop than the procedure outlined in the previous section at one-loop. Just as before, we will still fix the integrand on color-dressed cuts, however, now we can no longer simply decompose the full color dressed amplitude in terms of integrated color-stripped objects. At two-loop, a convenient basis of adjoint color factors happens to be the double-box and the cross-box, shown below: By iteratively applying Jacboi relations on internal edges, any color structure can be expressed completely in terms of functional relabelings of these two graphs. Concretely, these cubic graphs come dressed with the following color factors written in terms of group theory structure constants: The only functional difference between these color basis factors is the swap of the colored indices, b 1 and b 3 , in the second and third vertices. Just as at one-loop, we can carefully choose a spanning set of two-loop unitarity cuts in terms of our color-dressed NLSM operators in order to land directly on this minimal color basis. In doing so, we can focus our attention to integrating the kinematics that weight each of the color basis element. As we argued in the introduction, two-loop amplitudes for even-point theories, like NLSM, are completely determined by two recursive integrals, the double-bubble and the ostrichdiagram, modulo theory-dependent kinematic numerators. For NLSM, the relevant integrals are as follows: where the internal momenta appearing in the NLSM operators,p i = −p i andq i = −q i , are defined in terms of the two loop momenta as follows: Due to the simplicity of scalar theories like NLSM, the state sum for the double-bubble numerator is quite simple, and can be expressed concisely below: (4.26) where we have introduced the following notation above that combines internal and external momenta, τ (i) j = (k i + p j ) 2 and τ ij = (p i + p j ) 2 . Plugging in the color-dressed NLSM amplitudes, the same can be done for the ostrich-diagram cut, giving us .

(4.27)
The kinematic variables τ (i) j and τ ij are the same, except we have made the replacement p i → q i . The full two-loop NLSM amplitude can thus be computed by summing over the distinct labels of the resulting integrals, each weighted by suitable internal symmetry factors: We stress that the integrals appearing above are complicated by tensor integrals with many powers of loop momenta. However, as noted in the introduction, since both integrals are recursively one-loop, we can again apply a two-loop generalization of Passarino-Veltman for both the bubble and triangle integrals, stated in eq. (3.23) and eq. (3.38), respectively. With this, we proceed to the next step in the calculation.
Reduction Just as with one-loop, we will now reduce the integrals of eq. (4.28) to the a basis of D-dimensional scalar integrals using the EMU loop reduction. The double-bubble contribution yields the following: (4.29) Likewise, we can perform a similar reduction to the ostrich integrals -but with a small twist. The double-bubble integral above is completely separable, which allows one to treat the tensor reduction of each loop integral as independent. This is not so for the ostrich-diagram integral. The ostrich-diagram carries an internal bubble, I D 2 (q), with a loop dependent scale factor of q = l 2 + k 1 . Thus, we must first reduce the powers of l 1 and then perform a tensor reduction on all the remaining l 2 factors. Doing so yields the following result: (4.31) We have added the additional factor of (l 2 + k 1 ) 2 in the definition of [I 3 • I 2 ] D (k 12 ) to simplify the final expression in eq. (4.30). We will use this as our basis integral throughout the text. This two-loop integral can be evaluated analytically using the D-dimensional single scale bubble and triangle integral expressions in eq. (2.17). Thus, we have completely reduced the full two-loop NLSM amplitude to a linear combination of scalar integrals, I D 2 (k ij ) and [I 3 • I 2 ] D (k ij ). Now we will proceed to the final step of evaluating the amplitude in particular dimensions.
Exponentiation While there is plenty of physics to be extracted from the combination of eq. (4.29) and eq. (4.30), we will focus on the behavior of the leading divergence for D = 2−2ϵ for the target space model describing pions on CP 1 ∼ = SU (2)/U (1). The tree amplitudes for this theory can be obtained from eq. (4.4) by plugging in multi-trace color factors of the form, The delta function indices select out (anti)-holomorphic pion fields that live on the CP 1 target space. Plugging in explicit values for the "color" indices, we obtain the following tree-level amplitude: Using this tree-level contact as the seed in our EMU construction, we can directly compute the one-and two-loop amplitudes from eq. (4.12) and eq. (4.28). Furthermore, the leading divergences for the loop integrals in D = 2 − 2ϵ dimensions are the following: The analogous divergences for D = 4 − 2ϵ can be extracted directly from the integral expressions in eq. (2.17). Plugging these values in for the scattering process A(Z 1 ,Z 2 , Z 3 ,Z 4 ), we can extract the full leading logarithmic divergence from the one-loop color dressed amplitude in eq. (4.12), 36) and likewise for the two-loop result upon evaluating the integrals in D = 2 − 2ϵ, We write the two-loop amplitude in this suggestive form to emphasize that the logarithmic divergence present at one-loop in eq. (4.36) appears to exponentiate! More concretely, we have demonstrated through explicit calculation that through two-loop order, the leading divergence of the full color dressed amplitude goes as follows: This iterative form suggests that the IR divergence of the CP 1 model can be neatly extracted from the perturbative series as an overall factor, e Ω , where Ω is the divergent piece of the one-loop amplitude. Similar exponentiation of the IR has a long history in the context of gravity amplitudes [155], for which IR divergences appear due to soft graviton emissions that mediate long range interactions. Considering the established relationship [156][157][158][159] between coset manifolds (target spaces) and supergravity theories, there's a possibility that the exponentiation of the CP 1 model in eq. (4.38) is related to the analogous behavior of supergravity amplitudes [160][161][162][163][164].
In addition to the exponential structure of eq. (4.38) that we have computed for the CP 1 model, there is potentially a parallel story for the planar limit of the chiral NLSM amplitudes. Indeed, if one takes D → 2 in eq. (4.29), we find that the leading divergence of the double bubble in D = 2 − 2ϵ is precisely one half the square of the one-loop divergence in eq. (4.17). However, while the cross-box term appearing in eq. (4.30) is subleading when N c → 0, there are still relevant contributions from the double-boxes appearing in the ostrich integral. In principle one could add additional operators [165] or particle states [166] to eq. (2.22) that cancel the additional double-box contributions from the ostrich integral, while preserving the soft behavior needed for our EMU construction. Doing so would extend the exponential behavior to the subleading logarithms in the planar limit, similar to what was found in [167,168] for planar N = 4 sYM in D = 4 − 2ϵ. The iterated structure of planar N = 4 has been linked to the integrability of the theory [169][170][171][172][173][174][175][176][177][178]. We leave identifying such an extended theory consistent with exponentiation as a compelling direction of future study. Now that we have walked through our procedure for computing two-loop even point EFT amplitudes with NLSM as an exemplar, we are prepared to proceed to the main results of the text. While the construction of DBIVA integrands is more involved than constructing the simple expressions of NLSM, the general procedure is exactly the same. The only differences being additional gauge independent state sums, and the generally more complex tensor reduction.

DBIVA via EMU
In the remainder of this section, we will construct four-photon matrix elements in DBIVA theories through two-loop order. While we will consider a number of different internal matter states inside the loops, the interactions for which are not captured by the DBIVA Lagrangian eq. (2.35), we will require that all the external states are vectors. This will allow us to map our loop-level results to the basis of gauge invariant tensors described in section 3.3. Furthermore, in the last section, we will investigate how the loop level contributions turn on new operators in the EFT written in eq. (2.26).
As we noted in the introduction, DBIVA tree-level amplitudes, which appear in the field theory limit of the abelianized open superstring, can be realized as a double copy between NLSM and sYM [49], (4.39) To carry out this construction at the (multi-)loop level, we would simply replace the color factors of NLSM with color-dual loop-level numerators, and then integrate the result. In the case of N = 4 sYM, color-dual representations are known through four-loop [6,29]. Indeed, in section 4.3 we will use use double-copy to verify the N = 4 DBIVA results we construct here via unitarity. If we had access to color-dual representations for NLSM at two-loop, we additionally could plug those into the sYM integrands constructed from simplified methods of supersymmetric sums [179]. Just as in the previous section where we computed NLSM amplitudes, the recursively one-loop behavior will allow us to only consider the four-point contacts when constructing the physical parts of the two-loop integrand. We will start by reproducing the known one-loop results from a completely D-dimensional framework, and then move onto our novel two-loop results.

One-loop DBIVA
Analogous to our procedure in the previous section where we defined A NLSM (1|23|4) , first we will define a set of four-point operators for each distinct set of particle interaction. The Ddimensional contacts needed for DBIVA amplitudes at one-loop are as follows: where we have introduced notation (k a F b F c k a ) ≡ k µ a F µν b F νρ c k ρ a for the mixed scalar-vector amplitude. The particle content is labeled by γ for the BI photons, λ for the VA fermions, and X for the Dirac scalars. With these tree-level amplitudes in hand, we define the following set of matrix-elements needed for integrand construction: While these are a subset of the four-point operators we will use in the two-loop construction, they are sufficient for all the one-loop external-photon amplitudes. To construct the integrands, we will use the D-dimensional state sums for vectors and fermions, which we provide below: where q 2 = 0 is a null-reference momentum. Here, Γ µ are D-dimensional gamma matrices endowed with all the normal Clifford algebraic relations. Since we want to keep this as D-dimensional as possible, we define Γ 5 as the symbol representing the (D + 1)-th gamma matrix that anti-commutes with all other Γ µ . Using this definition of Γ 5 , we could also define a chiral projection state operator P ± = 1 2 (1 ± Γ 5 ) in a dimension agnostic form. Chiral operators become relevant for our external photon amplitudes when computing internal fermion contributions at two-loop. To see why, note that since we are only computing external photon amplitudes, the matrix elements must be parity even. This means that if there is a single chiral trace, the Γ 5 contribution must integrate to zero: (4.48) This property significantly simplifies the two-loop reduction in the presence of single chiral trace. However, at two-loop we can also get multi-trace contributions. In this case, the parity odd contribution (sourced by odd powers of Γ 5 ) must vanish after integration as follows: The second term is also parity even, and will contribute to D-dimensional Gramm determinants in the presence of two internal fermion loops. While in principle this term is relevant for two-loop diagrams with internal fermions, due to the simplicity of N = 4 DBIVA state sums, we won't need to account for it in our analysis at two-loop.
Before constructing the integrands, below we provide our conventions for the D-dimensional spinors and gamma matrices. We will normalize the gamma matrices as follows: (4.50) Furthermore, since we will eventually be evaluating our matrix elements in D = 4 after integration, we assume that the D-dimensional generalization of the Majorana condition holds throughout the calculation. That is, theū and v spinors obey the following relationship: where the D-dimensional charge conjugation operator can be defined in terms of the gamma matrices, Γ µ , as follows: From this, the spinor strings obey another in addition to the normal Clifford algebra identities [180]:ū This relationship essentially just imposes Fermi statistics on the identical Majorana fermions [181]. Now we are prepared to constructing the integrand with the state sums and operators defined above.
Construction Just as before, the first step in our procedure is to construct the integrand with all the internal loop dependence present. Taking the operators defined above, and applying the appropriate state sums, we obtain: 56) where exposed internal particles are taken to be on-shell. Above we have defined the internal momenta as q 1 = l and q 2 = −(l 1 + k 12 ). The vector and fermion contributions are rather complicated -however, to get a sense of what pops out of this D-dimensional construction, we provide the scalar cut below as an example: Constructing the full amplitude is then just a matter of summing over all cut contributions, each weighted by symmetry factors S α and the number, N α , of α-type particles: (2,3,4) . (4.58) Since the scalars are complex, and the fermions are oriented, they come dressed with symmetry factors of S X = −S λ = 1. The minus sign is due to the presence of a single fermion loop. Furthermore, since the photons are indistinguishable, they carry an internal symmetry factor of S γ = 2. Now we are prepared to carry out the integral reduction and evaluate each contribution to the one-loop amplitude above.
Reduction While the full D-dimensional integral at one-loop is rather complicated for the vectors and fermions, the tensor reduced integrals are actually quite simple -and its worth showing the result of applying eq. (3.23) to the integrands above in eq. (4.54). After reducing the tensor integrals that appear in eq. (4.58), we obtain (4.59) In terms of these integrals, the full amplitude at one-loop can be expressed as a sum over all distinct permutations of external legs, with each graph weighted by their associated symmetry factors: 4 (4.60) Now we can proceed to evaluating the integrals in the dimension of interest.
Integration While the procedure that we have employed thus far allows us to extract all of the logarithmic dependence of the loop integrals, for the remainder of this paper we are going for focus our attention to the leading order divergences in ϵ. The purpose of doing so is to identify the class of photon effective operators that we must add to the Born-Infeld action in order to cancel off anomalous rational terms in the S-matrix. We will leave future analysis of the novel DBIVA amplitudes in this paper as a direction of future study.
In order to capture all the tensor structures that survive after integration, including the evanescent contributions, we will map our results to the basis of tensor structures identified in section 3.3: We remind the reader that the photon structures given above are defined in eq. (3.41). At one-loop mass-dimension, there are only four tensor structures that contribute: Plugging in D = 4 − 2ϵ to the above evaluated one-loop integrals, the EFT expansion above forms a basis for all the algebraic (non-transcendental) parts of the one-loop integral. Explicitly, the Wilson coefficients above take on the follow values when we expand the integral in D = 4 − 2ϵ dimensions, We remind the reader that while there are four spanning D-dimensional four-photon structures, they have a non-trivial null space in four-dimensions where only three structures span the non-vanishing space. Since all of the above Wilson coefficients come dressed with a leading 1/ϵ divergence for pure photon amplitudes, this indicates the emergence of a divergent evanescent operator that will contribute at two-loop order. We will see the consequences of this when calculating the anomalous matrix elements for pure BI at the next loop order.
Projection While the expressions above capture the full behavior of the photon amplitudes, they obscure the 4D physics captured by spinor helicity variables. Rather than plugging in explicit values for D in the integral, it can be more informative to first project down onto a 4D basis of states, leaving the internal dimensional dependence untouched. Indeed, below we can see the effect of plugging in all-plus helicity configurations into the evaluated one-loop integrals above: (4.67) Here we can see two characteristic properties of the anomaly cancellation that takes place when we introduce scalars and fermions into our spectrum.
• First, all integrals carry a factor of (D − 4). This reflects the property that in D = 4, all the tree-level amplitudes that contribute to the cut must vanish outside the MHV sector. Thus, any contribution to the all plus matrix element must be suppressed by a factor of ϵ = (4 − D)/2 in order to push the logarithms to O(ϵ).
• Second, the common bracketed factor in all three integrals is weighted differently depending on whether the internal particle is a boson or a fermion. Bosonic contributions are weighted by the number of particles in the loop (in this case, 2 real scalars or (D − 2) helicity states), whereas the fermion weight is dependent on the dimension of the gamma matrix representation, where tr[Γ µ Γ µ ] = 2 D/2−1 D. In D = 4, all of these contribute equal magnitude to the full amplitude.
Putting this all together, we obtain the following non-vanishing algebraic parts of the one-loop matrix elements: The (− + ++) matrix element is identically zero due to D-dimensional four-point kinematics -the only available helicity structure is ⟨1|2|3] 2 [24] 2 which must be weighted by a massdimension 2 permutation invariant. The only such permutation invariant is s + t + u = 0, which vanishes regardless of the integration dimension. As a check, above we have reproduced the results found in [52], which used four-dimensional spinor helicity and dimensionshifting relations when performing the integration. This serves as a nice verification of our D-dimensional methods. Now we will proceed with the two-loop calculation.

Two-loop Born-Infeld
For all of the two-loop calculations, we will drop the header labelling which step in the process we are presenting. While we will omit these markers, our procedure is still the same: construction, reduction, integration and then projection. Due to the formidably large expressions that result from doing this calculation completely covariantly, most of our amplitudes will be presented after projecting down to 4D helicity states. To begin, we will compute the pure Born-Infeld two-loop amplitude. Like pions, there are only two diagrams that contribute at this loop order. Including symmetry factors, the full Born-Infeld amplitude at two-loop can be expressed as follows: As in the previous section, exposed legs implicitly sum over internal states. The grey blobs above represent D-dimensional t 8 F 4 operator insertions from the four-point BI tree amplitudes, which we labelled as M γγγγ (1234) . We give the internal loop momenta the same internal labels as we did the pions in eq. (4.26) and eq. (4.30), giving us the following integrals to evaluate for two-loop pure BI, (4.73) where the internal momenta obey the same convetion,p i = −p i andq i = −q i . Performing the tensor reduction on the internal loop momenta and projecting to 4D helicity states yields divergent quantities for all helicity configurations. In the interest of projecting to a Ddimensional basis of operators, and analyzing the U (1) anomaly present at two-loop, we will just focus on the leading divergences for each helicity configuration.
Leading (++++) divergence As we saw above, the one-loop matrix element has a rational all-plus contribution. This manifestly breaks the U (1) duality invariance present as tree-level. Moreover, this will contribute to non-vanishing 4D cut of the form: In principle, this divergence could be cancelled by the addition of a counterterm at O(α ′4 ) inserted into a one-loop matrix element with t 8 F 4 at O(α ′ 2 ). We will explore the effect of such anomaly cancelling counterterms in the next section, and we will see that this alone is not sufficient to cancel the divergence above. As we noted previously, there are additional 4D operators that appear at one-loop that will vanish when plugging in (+ + ++) physical states. These too could secretly contribute to the divergence expressed above. We can see this more clearly for the (− + ++) result below.
Leading (− + ++) divergence The same cut construction above suggests a different story for (− + ++) at two-loop. Indeed, the following cut should vanish when taken on-shell in D However, we found that M BI,1-loop (−+++) = 0 due to D-dimensional four-point kinematics. Despite this, we find similarly to the all-plus helicity configuration above in eq. (4.75), the one-minus matrix element also carries a leading order 1/ϵ divergence: At first glance, this appears to be a violation of the Optical Theorem. However, this divergence is sourced by one-loop evanescent operators that vanish in D = 4, but which carries a 1/ϵ divergence in D = 4 − 2ϵ. In the section 5.1, we will demonstrate this in more detail, and show how higher derivative four-photon operators can be used to cancel the U (1) anomaly.
In doing so, we can construct a D-dimensional quantum effective action that satisfies duality invariance in D = 4 through two-loop order in perturbation theory.
Leading (− − ++) divergence Finally, we express below the leading divergence for the aligned-helicity matrix elements. After reducing to the two-loop scalar integral basis, and projecting along 4D helicity states, we obtain the following expression at leading order in the ϵ-expansion:  [34] 2 helicity sectors are asymmetrically weighted in pure photon amplitudes. In the next section, we will show that these two operators carry equal weight when introducing additional states consistent with maximal supersymmetry.
Moreover, adding additional scalar and fermion states consistent with supersymmetry is the simplest way to cancel the U (1) anomaly computed above outside of the aligned helicity sectors. When summing over superfield contributions at loop level, the U (1) symmetry at treelevel is promoted to a U (1) R symmetry, and thus is protected pertubatively by supersymmetric Ward Identities. We demonstrate this explicitly in the next section by computing the two-loop four-photon matrix element in N = 4 DBIVA via our D-dimensional unitarity methods.

Two-loop N = 4 DBIVA
Much of the complication in performing a two-loop calculation in pure BI theory is the proliferation of high power of loop momenta left over after the state sum. These factors of loop momenta not only conspire with external momenta, but also mingle with external polarizations. There are a number of integral reduction algorithms [132][133][134][135][136][137][138][139] that are very effective for reducing factors of (k i · l), since they can trivially be expressed as a linear combinations of inverse propagators, However, the factors of (ε i ·l j ) can be more tedious, as they do not permit an inverse propagator expansion. This is part of the motivation for applying Passarino-Veltman to the recursively one-loop integrals present in the two-loop Born-Infeld amplitude. Luckily, the state-sewing for maximal supersymmetry is dramatically simpler than less than maximal supersymmetry [3,182]. This is most easily seen by considering the covariant operators we defined for one-loop DBIVA amplitudes. When one applies the conditions for maximal supersymmetry, stated below, then the state sum is completely independent of loop momenta and precisely reproduces the supersymmetric matrix element s 2 12 (t 8 F 4 ) when cut along the s 12 -channel. This statement holds covariantly in general dimension. Concretely, we find that for theories with maximal supersymmetry: This can similarly be used in the integrand construction N = 4 super-Yang-Mills, which has been computed to six-loop order [183]. In the case of maximal supersymmetry in D = 4, we can write the (t 8 F 4 ) (max) (1234) operator as a supersymmetric delta function with the all-plus permutation invariant introduced previously: where the delta function is a Grassmann valued polynomial of spinor-helicity variables: By applying this supersymmetric state sum, we find that the integrands for two-loop N = 4 DBIVA is trivially easy to construct -even simpler than NLSM integrands. Below we represent an internal on-shell superfield with a green line, and obtain the following integrand for the double-bubble: and similarly so for the ostrich-diagram integral contribution: (4.88) We should point out an additional verification can occur for any SUSY amount of DBIVA. We can recycle multiloop sYM integrands to verify our cuts. In general this is because double-copy works at tree-level, and we can always write the Jacobi-satisfying Yang-Mills cut numerators in terms of color-ordered sYM amplitudes, A sYM , offering us a KLT between cut-integrands. When we cut to four-point amplitudes, as we do for the cuts above, this is particularly simple as only one order needs contribute for each four-point tree, yielding:  Since there are no loop momenta interfering with the tensor structures of the external photons, the result of integration will just be proportional to a (t 8 F 4 ) tensor, which is manifestly U (1) duality invariant. Thus, as promised, adding in N = 4 states cancels the U (1) anomaly present in the pure Born-Infeld S-matrix. After applying the tensor reduction to the ostrich-diagram, we obtain the following expression for the two diagrams: This is easily evaluated in D = 4 − 2ϵ, from which we find the following expression for the leading order divergence of the four-photon two-loop amplitude in N = 4 DBIVA theory: With this result in hand, we will now demonstrate that the same calculation performed via generalized unitarity can be reproduced using loop-level double copy construction.

DBIVA via double copy
The calculation above was completely agnostic to (but verified by) the known tree-level relationship of DBIVA amplitudes as a double copy of NLSM and super Yang-Mills. As we will now show, the amplitude can be equivalently produced using the two-loop color-dual numerators of N = 4 sYM, and the NLSM integrands constructed earlier in the section. This observation provides further evidence for the consistency of double-copy construction at multi-loop order in perturbation theory, and serves as an existence proof for color-dual NLSM numerators at two-loop. At loop-level, the double-copy construction amounts to replacing the color factors in a cubic-graph representation of the integrand with a set of color-dual numerators. That is, starting with an L-loop NLSM integrand of the form: we can construct N = 4 DBIVA by replacing the color factors, C g , with the kinematic numerators, N N =4 g , of sYM amplitudes: where the kinematic numerators depend on particular choice of generalized gauge. In order for this construction to work, the kinematic numerators on at least one side of the double copy must satisfy all the same algebraic relations as the color factors. Since generalized unitarity allows us to construct the integrands on-shell, this construction conjecturally holds at loop-level as long as the gauge-theory is tree-level color-dual. While there are currently no color-dual NLSM numerators identified in the literature beyond one-loop, there are color-dual representation available for N = 4 sYM through fourloop four-point [31]. Below are the one-and two-loop basis numerators relevant for our construction: 99) and similarly so at two-loop four-point for the double-box and cross-box, which take on identical expressions: These two-loop numerators can be constructed via the rung-rule [182] and the no-triangle hypothesis [184]. We can see that both the one-and two-loop kinematic numerators are completely independent of loop momenta. This allows us to take our NLSM amplitudes from earlier in the section, and simply replace the color factors post integration. We will show the results of this procedure for both one-loop and two-loop DBIVA in the proceeding sections.

One-loop N = 4 DBIVA
At one-loop, performing the double copy is a simple task. As we noted above, since the box numerator for N = 4 sYM has no kinematic dependence on internal loop momenta, we can pull it out of the integration process. Thus, making the following replacement This is in agreement with the leading divergence of N = 4 DBIVA previously found in [52], along with our result in eq. (4.70) when taking N γ = 1, N λ = 4 and N X = 3, in concordance with D = 4 maximal supersymmetry

Two-loop N = 4 DBIVA
At two-loop it sufficient to show that replacing the color factors reproduces the integrals we found above with the unitarity construction. Starting with the NLSM integrands in eq. (4.29) and eq. (4.30), we make the same replacement of color factors post integration as was done at one-loop. We carry this out first for the double-bubble, and find that the full D-dimensional NLSM integral simplifies dramatically: (4.104) The cancellation between dimension dependent factors is even more startling for the ostrichdiagram integral: (4.105) This alone is sufficient to demonstrate the validity of the double-copy, as these integrated quantities are exactly the same as those produced via generalized unitarity of N = 4 DBIVA at two-loop in eq. (4.87) and eq. (4.88). We emphasize that based on our analysis, the doublecopy at two loop clearly holds in any dimension, as all the D-dependent prefactors drop out when replacing the color factors with N = 4 numerators. Considering the consistency of these two construction, this also serves as strong evidence for the existence of color-dual representation for two-loop pion integrands. We leave identifying such valid representation as an enticing future direction worth investigation.

Effective Actions
In this section we demonstrate how our basis of higher-derivative four-photon operators in eq. (3.41) can be used to construct quantum effective actions that capture loop-level effects. We will distinguish between position space operators that appear in the Lagrangian, O, and their corresponding matrix elements, T , which are on-shell quantities in momentum space: Similar to the matrix elements, the operators can be expressed in both D-dimensional and 4D representations. For example, we can define O F 2 F 2 operators as follows: where the spacetimes indices run from 1, 2, ..., D, and F µν = ∂ µ A ν − ∂ ν A µ are operator valued abelian field strengths. These operator valued expressions can be normalized so that for an all outgoing four-photon scattering process, |out⟩ = |k 1 , k 2 , k 3 , k 4 ⟩, they are related to the on-shell tensor basis elements defined previously in eq. (3.41) as follows: Likewise we can construct field theory operators for the split helicity operators using the spinor form of the field strengths: where the dual field strength are defined asF µν = 1 2 ϵ µνρσ F µν , and the spin matrices are expressed as σ µν ∓ in terms of the Pauli matrices as follows: In the above expressions, the spacetime indices now run from 1 to 4. Using these manifestly 4D operators, we can similarly define operators that select out particular helicity states and derivative structures: where the trace is over the spinor indices of the Pauli matrices. As defined, there operators will produce the 4D helicity structures defined in eq. (3.43) and eq. (3.44) when contracted between the outgoing helicity state, |out⟩ = |k − 1 , k − 2 , k + 3 , k + 4 ⟩, and the incoming vacuum state: In the duality invariant quantum actions we construct in this section, we will reserve O for effective actions, and T for on-shell matrix elements used in the EMU construction. However, since the field theory operators are redundant up to field redefinition and equations of motion, we'll find it convenient to work in terms of the on-shell matrix elements, and then implicitly define the operators in terms of these on-shell expressions.
Below we first study anomaly cancellation for the multi-loop photon amplitudes computed above. After this, we proceed by interpreting these higher derivative operators as double copies between NLSM and higher derivative Yang-Mills amplitudes with off-shell higher-spin modes. As was demonstrated in previous work by the authors [51], the full set of D-dimensional four-photon operators can be constructed via adjoint double-copy, but at the cost of introducing off-shell higher spin modes in the single-copy vector theory. However, these higher-spin modes can be absorbed consistently in a double-copy framework by introducing symmetric algebraic structures. We discuss future applications of this construction in section 6.

Anomaly cancellation
Now we begin with our study of higher-derivative extensions to BI theory. Our goal is to identify the higher derivative four-photon operators that are needed to cancel the U (1) anomalous matrix elements computed in the previous section. We will start with the one-loop corrections, captured by tree level insertions at O(α ′4 ), and then proceed with the two loop corrections, which combine both O(α ′4 ) operator insertions at one-loop and O(α ′6 ) operator insertions at tree-level. In doing so, we demonstrate that cancelling the U (1) anomaly through two loop order can be achieved with local finite counterterms if and only if we introduce an evanescent operator at O(α ′4 ) to the Born-Infeld action.

One-loop
In section 4.2 we computed the one-loop matrix element for a general DBIVA theory. Plugging in the values N γ = 1 and N λ = N X = 0 we obtain the following anomalous all-plus matrix element for pure Born-Infeld theory: In Ref. [52], Elvang, Hadjiantonis, Jones, and Paranjape identify a candidate 4D counterterm that cancels this anomalous matrix element, whose prediction we have called T 4D (++++) , thereby restoring duality invariance through one-loop. As noted in the previous section, the one-loop matrix element can be mapped to our D-dimensional operator basis. In general, all available 4D tensor structures at O(α ′4 ) map onto our D-dimensional basis. One particular map we provide below: where we have defined the following D-dimensional operator that projects down to the all-plus configuration, and all the 4D operators have the freedom to add the previously defined evanescent operator, T ev. , Thus, we can construct the new effective photon Lagrangian, L BI+CT , with the addition of the all-plus counter-terms to our Born-Infeld Lagrangian: where we have used eq. (5.1) to implicitly define the operators appearing in the quantum effective action above. While there are a number of perturbatively equivalent construction of these operators, below we provide a couple expressions that resemble the on-shell basis elements: Thus, eq. (5.18) constitutes a duality invariant photon theory through one-loop order. With this, we can identify what additional operators will be needed to cancel the anomaly through two-loop.

Two-loop
The first step in identifying the requisite operators needed to cancel the two-loop anomaly at O(α ′6 ) is to perform another one-loop calculation at this mass dimension, which includes the counterterms of L BI+CT defined above. The one-loop amplitude is constructed as follows: Both of these operator insertions can be evaluated using the same D-dimensional procedure used throughout the text. The all-plus counterterm yields the following contributions to (+ + ++) and (− + ++) helicity configurations: We note that there is a distinction between the first and second (−+++) expressions. The first expression is dressed with (D − 4) 2 , which pushes the leading contribution to O(ϵ). Whereas the second term is identically zero because the 4D helicity structure carries an overall factor of (s + t + u) = 0. In addition, since there is a non-vanishing 4D residue for the all-plus integrand, the integral has a leading order divergence in ϵ.
Below we find it instructive to show the D-dependence of the evanescent operator insertion, which yields the following matrix element contributions: This will produce O(ϵ 0 ) matrix elements in the (− + ++) helicity sector. Thus, in order to cancel the divergent part of the two-loop (− + ++) anomaly computed in eq. (4.77), we must weight the evanescent operator by a numerical factor that diverges in D = 4. Given the particular numerical value computed in the previous section, we find the evanescent Wilson coefficient must take the following D-dependent value: where the (D − 4) in the denominator cancels the factor in the numerator above. In order to further absorb the remaining rational terms, we must introduce an additional set of tree-level operators at O(α ′6 ). At this order in mass-dimension, there are seven distinct operators: By adding these operators to the effective Lagrangian above, we have verified that there is sufficient freedom to absorb the remaining rational terms present at two-loop. Of these available operators, only the F 3 tensor structure is non-vanishing when projected along the (− + ++) helicity configuration. Furthermore, just as at O(α ′4 ), there is a single evanescent matrix element which we define below: Thus, of the seven available D-dimensional operators, they are projected to only six distinct 4D tensor structures. We will describe the counting of 4D versus general dimension photon operators in generality at all orders in α ′ in more depth at the end of this section.

Double copy construction
As we have stated in the text, it is well known that DBIVA theory can be constructed at tree-level as an adjoint double copy between NLSM and sYM amplitudes. There is now a large body of literature studying double-copy construction of higher derivative gauge theory counterterms [13,[15][16][17][18][19][21][22][23][24], like those used above to cancel U (1) anomalous matrix elements in pure Born-Infeld theory. Indeed, recent work by the authors demonstrated that all four-photon operators can be constructed consistently via the double-copy [51]. Here we briefly describe the single-copy gauge theory that when double copied with NLSM produces the higher derivative operators of the previous section.

Symmetric-structure double-copy
To realize the double-copy construction that produces the counterterms above, Ref. [51] first decomposed NLSM pions amplitude into symmetric structure constants using the U (N ) color identity where the symmetric structure constant is defined as follows By applying this color algebra identity to the four-point NLSM amplitudes of eq. (4.4), one finds that pions can similarly be expressed as a symmetric-structure double copy: where c dd s ≡ d abe d ecd and the NLSM symmetric s-channel numerator is n dd,π g = s 2 . By identifying a set of gauge theory numerators that obey the same algebraic relations as the color factors, one can construct consistent double copy theories.
For example, consider the two-loop divergence for the (++++) anomalous matrix element in pure BI theory. The 4D helicity structure can be captured by a symmetric structure double-copy between NLSM pion numerators and a local gauge theory contact at O(α ′5 ). The s-channel numerators for this symmetric-structure double copy are as follows While this construction lacks any algebraic relations between the four-point kinematic-factors, similar symmetric numerators were found at six-point for NLSM, which obey non-trivial algebraic relations [51]. In addition to the photon counterterms constructed above, symmetric double copy is likewise a natural description of gravitational counterterms. Rather than composing the local vector numerators, n dd,HD , with symmetric NLSM numerators, n dd,π g , they can also be squared -yielding a gravitational counterterm: With these building blocks, symmetric double copy construction captures the exceptional four graviton amplitude of Ref. [185], which only considered linear combinations of gauge theory amplitudes rather than gauge theory numerators. We note that a similar construction was recently identified in Ref. [186] beyond four-point, which demonstrated that gravity amplitudes permit double copy construction in terms of gauge-invariant numerators that do not obey any algebraic relations between themselves, much like the symmetric numerators above in eq. (5.35). It would be fascinating to determine if the resulting adjoint higher-spin theory that is sourced by the anomolous matrix elements in pure BI theory is in any way a consistent physical theory. One natural guess would be some sort of tensionless limit of string theory, where the tower of massless higher spin modes are needed for the theory to be unitary. Along these lines, there have been many recent studies into building higher-spin gauge theories constructively from general principles of locality and unitarity [187][188][189][190][191][192]. There is also a possibility that these higher-derivative operators could be required for double-copy consistency at highermultiplicity. Such behavior was demonstrated for YM + F 3 gauge theory [18], where doublecopy consistency demanded a tower of four-point contacts that conjecturally resums to DF 2 + YM theory of Johansson and Nohle [193,194]. More recently, similar structure was likewise identified for higher-derivative corrections to biadjoint scalar EFT [22,23]. We see better understanding the double-copy consistency of higher-spin gauge theory as an exciting direction of future study.

Evanescent operator counting
Before concluding, in this section we compute the counting of evanescent operators at higher orders in mass dimension needed for U (1) duality satisfying quantum effective actions. We will do so using Hilbert series, which are an effective method for counting the number of independent operators at a particular order in mass dimension [195,196] and have been used in abundance in the SMEFT literature [197][198][199]. Hilbert series account for Lagrangian-level redundancy by counting the number of distinct on-shell quantities. For example, the following two operators would only be counted once since they are equivalent up to a boundary term, due to the Bianchi identity, D (µ F νρ) = 0. As such, we will use the on-shell matrix elements, T , for the computations in this section. Below we construct the Hilbert series for general dimension photon operators, H gen.D , and similarly so for D = 4 operators, H D=4 , for which we are restricted to a subset of 4D helicity operators. To construct the Hilbert series for photon operators, we will define the coefficient at order-n of a polynomial in α to be the number of operators that appear at O(α ′n+2 ) in our dimensionful coupling. It turns out that there are only two integer sequences that control the operator counting, which are common feature when enumerating operator bases for 2-to-2 scattering events [200][201][202][203][204][205]. The first sequence is produced by Hilbert series that counts four-point permutation invariants, σ x 3 σ y 2 , which we call H (ijkl) , [H (ijkl) ] = 1, 0, 1, 1, 1, 1, 2, 1, 2, 2, 2, 2, 3, 2, 3, 3, 3, ... H ev. = α 2 (α − 1) 2 (α 2 + α + 1) . (5.53) It would be interesting to determine whether the three Hilbert series stated above owe their construction to some hidden geometric origin. Indeed all of the operator counting in the SMEFT literature can be traced to the geometry of the group theory representations that underly the Standard Model [195,196]. We leave identifying these concealed mathematical structures as an exciting direction of future study.

Conclusions
In this manuscript, we have carried out a detailed analysis of even-point effective field theories through two-loop in the perturbative expansion. In section 2, we provided a review of the generalized unitarity and integration methods employed throughout the text. Then in section 3 we introduced and developed the on-shell constructive method of Even-point Multiloop Unitarity (EMU), and computed the tensor reduction for triangle and bubble integrals in D-dimensions of arbitrary rank along with a spanning set of four-photon operators needed to capture higher-loop effects. Due to the simplicity of even-point multi-loop amplitudes of the nonlinear sigma model (NLSM) and Dirac-Born-Infeld-Volkov-Akulov (DBIVA) theories, these methods allowed us to compute fully integrated two-loop amplitudes for NLSM, pure Born-Infeld, and N = 4 DBIVA theory in section 4. Finally, in section 5 we studied the quantum effective actions that capture the aforementioned loop effects, and studied the all-order counting of higher-derivative four-photon operators using Hilbert series. In doing so, we have identified a variety of rich physical structures that we summarize below: Exponentiation In eq. (4.29) and eq. (4.30) we computed in general dimension, D, the two contributions to NLSM two-loop amplitudes. Plugging in explicit color structures, we found that the leading divergence of NLSM amplitudes on a CP 1 target space exponentiate in D = 2 − 2ϵ, akin to the IR exponentiation of gravity amplitudes [155,[160][161][162][163][164]. In addition, evaluating these diagrams in the planar limit, where N c → ∞, we found that the doublebubble integral iterates the leading divergence and subleading logarithms in the full one-loop amplitude. Thus, identifying corrections to the NLSM Lagrangian that absorb the ostrich diagram would imbue the scale-dependent logarithms with exponential structure at loop-level. This non-trivial property is found in theories that are conjectured to be integrable [206][207][208], like N = 4 super-Yang-Mills [167,168] in D = 4 − 2ϵ. In future work, we hope to study whether this iterative structure can be further applied beyond the CP 1 model.
Anomalies Equipped with our D-dimensional integration methods, we performed a similar calculation at two-loop for pure Born-Infeld in section 4.2.2. In doing so, we demonstrated that the previously identified one-loop counterterm [52] is not sufficient to cancel the two-loop anomaly. In fact, the (− + ++) anomaly of eq. (4.77), which was absent at one-loop, diverges at two-loop order due to the presence of a one-loop evanescent operator in the D-dimensional formulation of Born-Infeld theory. One resolution to this anomaly comes in the form of introducing N = 4 DBIVA superfields in the two-loop state-sum. This protects the S-matrix from anomalies by promoting the classically conserved U (1) duality to a supersymmetric Rsymmetry. The evaluated integrals that contribute to the N = 4 DBIVA two-loop amplitudes are provided in eq. (4.87) and eq. (4.88).
As mentioned at the close of section 2.5, higher-multiplicity tree-level abelianized opensuperstring amplitudes, in four-dimensions, violates U (1) duality at higher-derivative corrections beyond the leading DBIVA predictions. It would be fascinating if these higher derivative violations of U (1) duality in the OSS spectrum are precisely those needed to cancel anomalous behavior at the quantum level. To test this would require computing a six-point one-loop amplitude in supersymmetric DBIVA. Such a calculation could in principle be performed with the double copy using the color-dual integrands recently constructed in Ref. [34]. We see this as a natural future direction and application of the methods we have developed here.
Evanescence Another resolution to the two-loop anomaly comes in the form of higherderivative pure-photon counterterms. As we demonstrated in section 5.1, one must introduce a divergent evanescent operator at one-loop in order to absorb the two-loop anomaly. This is similar to the anomalies of pure Einstein-Hilbert gravity [145,146], which vanish at oneloop since the R 2 Gauss-Bonet term is evanescent in D = 4, but which diverge at twoloop order [60][61][62]. We have constructed the higher derivative Lagrangian through O(α ′6 ) in eq. (5.29) that cancels the divergent part of the anomaly, along with a spanning set of counterterms in eq. (5.30) needed to absorb the finite part left over at two-loop. Given the complexity of gravity calculations at high loop order, further studies of multi-loop Born-Infeld amplitudes could serve as an accessible laboratory for studying evanescent effects beyond one-loop in double-copy constructible theories. To this end, we have used Hilbert series to count the number of four-photon evanescent operators to higher-order derivative corrections in section 5.3 to aid in future studies.
In addition to these themes woven throughout the text, we have, en passant, identified novel double-copy structures at two-loop. In section 4.3 we found that two-loop N = 4 DBIVA amplitudes can be constructed via the double copy of color-dual N = 4 sYM integrands with the generalized unitarity cuts of NLSM. This construction was D-dimensionally identical to the result obtained via generalized unitarity and maximal supersymmetric state sums. While not a proof, this provides strong evidence for the compatibility of NLSM with color-kinematics duality beyond one-loop. This result also serves as the first non-gravitational double-copy beyond one-loop, further supporting the consistency of color-kinematics duality at loop-level, which as of today remains a conjecture.
Moreover, recent work by the authors has demonstrated that color-kinematics duality can be used as a bootstrap principle to constrain higher derivative operators. This has been shown both for gauged NLSM amplitudes [19] and YM + F 3 theory [18], the later of which is particularly relevant for anomaly cancellation, and possibly UV completion, of both N = 4 supergravity and the R 3 modification to Einstein-Hilbert gravity [18]. Indeed, similar structure has been recently identified in color-dual scalar theories [22][23][24]. This observation suggests a new paradigm that elevates color-kinematics duality from a mathematical correspondence capable of encoding IR symmetries, to a principle that probes signatures of UV physics captured by higher-derivative corrections. Guided by this new paradigm, a natural next step is to determine whether the anomaly cancelling counterterms of eq. (4.75) and eq. (4.77) source additional higher-loop counterterms constrained by double-copy consistency, in the spirit of [18]. We see this as an exciting future direction in further understanding the loop-level constraints imposed by the duality between color and kinematics.