Rigorous Bounds on Light-by-Light Scattering

We bound EFT coefficients appearing in $2 \to 2$ photon scattering amplitudes in four dimensions. After reviewing unitarity and positivity conditions in this context, we use dispersion relations and crossing symmetry to compute sum rules and null constraints. This allows us to derive new rigorous bounds on operators with four, six, and eight derivatives, including two-sided bounds on their ratios. Comparing with a number of partial UV completions, we find that some of our bounds are saturated by the amplitudes that arise from integrating out a massive scalar or axion, while others suggest the existence of unknown amplitudes.


Introduction
Many physical systems display multiple characteristic energy scales, where they are described by different degrees of freedom and governed by different dynamics. When these scales are well separated, it is possible to integrate out the high-energy excitations and restrict attention to low-energy ones, which will be described by a new set of interactions. The resulting Effective Field Theory (EFT) represents a powerful framework to parametrize the ignorance about the microscopic behavior of a model. Instead of worrying about the detailed ultraviolet (UV) completion of a theory, one can instead specify the low-energy content of the theory and expand the interaction in a controlled series of higher-dimension operators made with a few building blocks and suppressed by increasing powers of the microscopic scale. The coefficients of this expansion, also called Wilson coefficients, encode the UV details of the system. It was realized some time ago [1][2][3] that not all values of these coefficients are consistent with a well behaved UV completion: despite the spirit of the EFT approach to be as agnostic as possible about the microscopic description of the theory, there are essential properties we should not abandon if we want the underlying theory to be consistent with basic axioms obeyed by S-matrices. Enforcing these conditions on the EFT produces a set of linear inequalities involving the Wilson coefficients.
principle, as the constraints are very simple when expressed in terms of EFT coefficients. The two-sided bounds obtained for EFT coefficients are reminiscent of the islands of allowed values for operator scaling dimensions and couplings in the conformal bootstrap. As we will discuss, the analogy extends even further: we find examples of known UV completions appearing at kinks on the boundary of the allowed region. This raises the exciting possibility that the EFT bootstrap could be an efficient method for finding unknown UV completions for a given set of low-energy fields.

Photons EFT and summary of results
In this paper we continue the exploration of EFT constraints by focusing on photons in four dimensions, i.e. massless spin one vectors. Neglecting gravity, the photon is the only massless field in the Standard Model of particle physics which is not confined at low energies. Hence the electromagnetic gauge field A µ is the only propagating degree of freedom in the infrared (IR). An observer performing experiments at center of mass (COM) energies √ s much lower than mass of the lightest (charged) particle will only be able to scatter photons and the fundamental observables are therefore the scattering amplitudes A λ 1 λ 2 ···λn , where λ i = ± denotes the polarization of the i−th photon. Such an observer might try to describe the outcome of the scattering experiments by an EFT with a Lagrangian of the form where a 1 , a 2 are dimensionful constants (Wilson coefficients) andF µν = 1 2 µνρσ F ρσ . Dimensional analysis implies that the importance of the Wilson coefficients grows with COM energy, meaning that our low-energy observer will conclude that "new physics" must appear at some energy scale M . We will refer to this new physics as a partial UV completion, to emphasize that it may itself need further completion at some even higher energy scale.
In the Standard Model, the scale M is determined by the electron mass M 2 = 4m 2 e , and the leading correction to the Wilson coefficients come from integrating out an electron loop, giving (1.2) With these values, (1.1) agrees with the famous Euler-Heisenberg Lagrangian, written down in the 1930's as a precursor to quantum electrodynamics (QED) [15][16][17]. 1 At higher energies, the contributions from heavier charged particles of the Standard Model become important [20][21][22]. Additional contributions coming from physics beyond the Standard Model, such as a light axion, could alter these predictions, see section 4.4.1. This has triggered several experimental efforts to measure deviations from (1.2). The Wilson coefficients in the Euler-Heisenberg Lagrangian may be observed in a magnetic birefringence experiment, where a strong applied magnetic field gives rise to an anisotropy in the optical index proportional to a 2 − a 1 . 2 In addition, in a direct scattering experiment γγ → γγ of unpolarized beams, the total cross-section at leading order is proportional to a 2 1 + a 2 2 − 2/3a 1 a 2 (see e.g. [24]). More generally, we can ask what partial UV completions are allowed based on some very general assumptions such as unitarity, causality and analyticity of scattering amplitudes. This question has a long history in the case of a low-energy EFT of scalars [1][2][3][25][26][27][28]. The case of photons was investigated in [29], and massive vectors in [30]. 3 These papers were able to derive positivity conditions on various Wilson coefficients, "one-sided bounds". In particular, [29] showed a 1 + a 2 0 using similar arguments involving dispersion relations in the forward limit, and a 1 0, a 2 0 individually using independent arguments involving the absence of tachyons and ghosts plus some stronger assumptions about the form of the UV completion. Later efforts [34][35][36] showed that both a 1 and a 2 are individually positive based on the dispersion relation arguments of [29] alone. In this paper, we use the method of [7,37] to show that it is possible to get two-sided bounds on the Wilson coefficients for low-energy electromagnetism. In particular, we show how crossing symmetry leads to redundancy in the low-energy description, allowing the derivation of null constraints, which are non-trivial constraints on the high-energy partial wave densities. These drastically strengthen the possible bounds, especially when combined with semidefinite programming.
To give an idea of the various bounds we find, consider the following parametrizations of the low-energy amplitudes: 4 where the constants g 2 and f 2 are related to Euler-Heisenberg coefficients simply by In the present work, we find a number of further interesting bounds, including: • A two-sided bound for the ratio g 3 /g 2 of the form −4.828427 < g 3 M 2 g 2

1.
The upper bound is optimal, while the lower bound is weakly sensitive to the numerical precision (number of null constraints), to be specified below.
• A finite allowed region in the space spanned by the parameters figure 4. In this parameter space, we compare our bound on (g 4,1 , g 4,2 ) to the one found in [4].
• A finite allowed region in the space of couplings f 2 , f 3 , f 4 , g 3 , g 4,1 , g 4,2 , normalized to g 2 and multiplied by the appropriate power of M . Sections of this region are shown in figures 2, 3, 4, 5 and 6.
Perhaps more interestingly, we find that a number of our bounds are saturated by amplitudes that arise from integrating a single massive scalar or axion at tree-level. We find that this is not the case for amplitudes arising from integrating out a massive graviton. Thus the scalar and axion tree-level exchange amplitudes are extremal, and those of graviton exchange are not. The outline of the paper is the following: In section 2, we show explicitly how to parameterize the low-energy (EFT) amplitudes in terms of symmetric polynomials of Mandelstam invariants, and the high-energy amplitudes in terms of partial waves with unknown but positive spectral densities. In section 3, we derive sum rules and null constraints by applying dispersion relations to these amplitudes, and use use them in section 4 to derive both analytic and numerical bounds on the EFT coefficients. Section 4 also contains a comparison with partial UV completions derived from integrating out massive particles at tree-level or one-loop level. We finish with a discussion and some appendices providing more details on some aspects mentioned in the main text.

Set-Up
The method we use in this paper can be summarized in three steps: 1. parametrize the low-energy amplitude as EFT coefficients times Mandelstam invariants, 2. parametrize the high-energy amplitude by partial waves, 3. use dispersion relations to relate the low-energy and high-energy amplitudes.
In this section, we explain how we perform the first two parts of this strategy.

Four-photon amplitudes
We would like to understand: what effective field theories of photons are consistent with unitarity, Lorentz invariant S-matrices in the ultraviolet? To begin to answer this question, we will consider scattering amplitudes with four external photons. Such amplitudes may be written using a basis of 16 different observables, A λ 1 λ 2 λ 3 λ 4 , where λ i = ± is the helicity of particle i, and we use all ingoing conventions. These amplitudes are related by parity, time-reversal, and boson exchange according to: In what follows, we will assume that our theory satisfies all three of these symmetries. Together, P, T , and B reduce the original basis of 16 amplitudes to only five independent ones. We may further consider crossing symmetry, which leads to the following requirements: A ++++ (s, t, u) : fully s-t-u symmetric, Crossing symmetry therefore allows us to write the amplitudes in terms of only three functions: (2.5) This will be the low-energy expansion of our theory. The functions f , g, and h are polynomials of the Mandelstams given by f (s, t, u) = f 2 (s 2 + t 2 + u 2 ) + f 3 stu + f 4 (s 2 + t 2 + u 2 ) 2 + . . . , g(s|t, u) = g 2 s 2 + g 3 s 3 + g 4,1 s 4 + g 4,2 s 2 (s 2 + t 2 + u 2 ) + . . . , h(s, t, u) = h 3 stu + h 5 stu(s 2 + t 2 + u 2 ) + . . . .

(2.6)
These functions are fixed by the symmetries. f includes every term which is totally symmetric in s, t, and u. h is also symmetric, and includes every monomial which is stu times a symmetric term. g(s|t, u) is s 2 times every monomial which is symmetric in t and u. They represent non-renormalizable interactions characteristic of the EFT picture we use here. The Lagrangian that leads to a particular amplitude is not unique in general because field redefinitions may change the Lagrangian but may not change the amplitudes. As an example, the amplitudes (2.6) can be found from quoted in the introduction. In general, the constants f k , g k(,i) and h k will multiply operators containing 2k derivatives, and we will therefore refer to these constants as 2k-derivative coefficients.

Partial wave expansion
Next in this analysis is to understand which of these helicity amplitudes can be used to derive constraints. This hinges on them having a partial wave expansion with certain positivity properties. We may, of course, take linear combinations of amplitudes -in fact, we will need to do so to ensure that our observables retain the s ↔ t crossing symmetry that is necessary for the contour deformation argument in the next section. The general method of analyzing unitarity constraints in the context of spinning external operators was nicely outlined in [38] (see also [34,39] for previous discussion on constraints from unitarity and crossing for external particles with spin). The partial wave expansion for spinning external particles with helicities λ i takes the form [40] where cos θ = 1 + 2u s , and where the Wigner d-matrices d a,b generalize the Legendre polynomials to spinning external states. The index i depends on the external helicities. Explicitly, the expansion for each of our five amplitudes reads:

Sum Rules from Dispersion Relations
We have now established the low-and high-energy parameterizations of our amplitudes using the EFT and the partial wave expansions, respectively. Dispersion relations relate these by requiring that a contour integral over the entire amplitude in the complex s-plane vanishes. This requires the following assumptions: 1. Causality: for fixed u < 0, A(s, u) is analytic on the upper half-plane, Im(s) > 0.

2.
Regge boundedness: for fixed u < 0, the total amplitude falls off faster than s 2 for large s. Specifically, This large-s behavior has been established rigorously for scattering in gapped theories [42,43], however in our case we will take it as an unproven assumption.
3. Weak coupling: we assume that the low-energy amplitudes are weakly coupled at least up to the scale of new physics M , so that low-energy loops are suppressed. As a result, M is the lowest energy where cuts could possibly appear in the amplitude. In general, we will get stronger bounds with a higher value of M , and weaker bounds with a lower, or "more conservative" value of M .
We define the amplitude on the lower half-plane by analytic continuation, i.e. A(s * , u * ) := A * (s, u). So the first assumption implies that the amplitude is analytic everywhere but the real s-line.
These assumptions imply what we call a doubly-subtracted dispersion relation: Doubly-subtracted refers to the two additional powers of s in the denominator compared to the factor (s − s ) −1 . Now we assume that, at low energies, the amplitude is given by our EFT expansion, and at high energies, it is given by the partial wave expansion. The integral in (3.2) obtains contributions from the poles at s = 0, s = s, and s = −u, and from cuts on the real-s axis, starting at s = M 2 and s = −M 2 − u. In this paper, we will consider only s − t symmetric amplitudes, which allows us to relate the s = s pole to the s = −u pole, and the left-handed cut to the right-handed cut (similar to the set-up in [37]). We visualize the analytic structure and the combination of cuts in figure 1. As a result, the dispersion relation becomes .  The contour is deformed inwards, picking up contributions from the three poles indicated and from two cuts on the real s axis, starting at M 2 and −M 2 − u respectively. Based on the assumed symmetry s ↔ t = −s − u, the integrals over the two cuts can be combined.

Sum rules
The fact that the contour integral isolates the imaginary part of the partial wave expansion is key to this argument. It means that we can derive positivity bounds on the low-energy amplitudes, and thus the EFT coefficients, by choosing observables which have positive partial wave expansions. We will consider general s ↔ t symmetric linear combinations of the amplitudes defined by which has a low-energy expansion The unitarity considerations in section 2 imply that A H [x 1 , x 2 ] satisfies appropriate positivity conditions for x 1 ∈ [−1, 1] and any x 2 . Let us see what sum rules result from considering these amplitudes.
Starting from (3.4) and using the dispersion relation (3.3) we get and where we have defined the brackets 1], all of these brackets represent an integral over a positive measure for all J in the respective sums. In the following, it will be convenient to define Now we are in a position to derive multiple sum rules from the "master sum rule" (3.6). We do this with a double expansion in the small s and small u limit. Due to the symmetry s ↔ −s − u, we follow [37] and replace the expansion in s with the more convenient expansion in s(s + u). At each order in s(s + u) we get a sum rule C 2n,u [x 1 , (3.10)

Four-, six-, and eight-derivative operators
To isolate the dependence on the low-energy EFT coefficients, we next expand the sum rules C 2n,u [x 1 , x 2 ] at small u, or equivalently take u-derivatives followed by setting u = 0. Let us look at the results of doing this for the first few orders in the s and u expansion: where J 2 = J(J + 1). Using these amplitudes, we can solve for the EFT coefficients in terms of the brackets. For example, In table 4 in appendix C.1, we collect the sum rules for all the seven EFT coefficients that are analyzed in this paper: f 2 , g 2 , f 3 , g 3 , f 4 , g 4,1 and g 4,2 .

Null constraints
The sum rules derived above contain redundancies, the first of which happens at eightderivative order. This allows us to derive null constraints by writing the same coefficient two different ways [37]. For example, which directly leads to the null constraint As we have removed f 4 from this equation, X 2,0 represents a constraint on high-energy data only. Enforcing this constraint allows for the derivation of stronger numerical bounds -essentially, instead of searching for functionals which a positive for all possible highenergy data, we can look for functionals which are positive for all possible high-energy data satisfying (3.15). We can do the same for the combination g 4,1 + 2g 4,2 , giving These null constraints ultimately derive from crossing symmetry, which restricts the number of possible EFT coefficients. 5 The result is that there are more sum rules than coefficients, hence the redundancy we see above. The null constraints inherently apply to only the high-energy part of the amplitude; in essence they bound the partial waves, and seem to imply a weak form of low spin dominance, similar to what is considered in [41].

Systematics for null constraints
Though we will only analyze the seven operators with eight or fewer derivatives, we will need to consider higher sum rules in order to generate more null constraints. In practice, it is possible to generate hundreds of null constraints from these formulas in Mathematica. Such null constraints are most straightforwardly computed by considering completely s-t-u symmetric amplitudes. We have two linearly independent amplitudes satisfying this property: and Crossing symmetry simplifies the low-energy expansions to We shall refer to null constraints that arise from the first amplitude (3.17) as "f -type", and those from the second amplitude (3.18) as "g-type". These both arise from completely st-u symmetric polynomials. It is easy to realize that the most general such amplitude can be written as a linear combination of simple polynomials (s 2 +t 2 +u 2 ) a (stu) b . The structure of the low-energy amplitude is therefore completely equivalent to the one considered in [7,37].
The high-energy expansion is different. In [37], general null constraints were written down on the form where P J (cos θ) are Gegenbauer polynomials (which reduce to Legendre polynomials in four dimensions). Since we have exactly the same cancellations on the low-energy side as in that paper, we can find f -type null constraints in the case of photon scattering by replacing Likewise, the g-type null constraints follow from Note that the P J (cos θ) = d J 0,0 (θ), so that the expressions entering in the 0 brackets are identical to those entering in the scalar bracket of [7,37].

Results
In this section, we report what comes out of the dispersion relation methods above. We are interested in analyzing all coefficients for operators with up to eight derivatives. We will display a variety of analytic and numerical bounds.

Bounds on four-derivative coefficients
Recall from section 3 that As a result, we find that Since all the involved brackets, defined by (3.8) and (3.9), are manifestly non-negative, we are left with two positivity conditions, or equivalently |f 2 | g 2 . This shows that the two constants a 1 and a 2 defined in (1.1) are individually non-negative: a 1 0, a 2 0. In [29], this conclusion was reached by assuming an ansatz for the UV completion and requiring the absence of tachyons and ghosts. The more general dispersion relation method employed in that paper only implied the result a 1 + a 2 0, equivalent to g 2 0. It was later observed using dispersion relations in [34,35].

Null constraints and bounds on six-derivative coefficients
Next we move to coefficients g 3 and f 3 appearing at six-derivative order. These are given by where J 2 = J(J + 1). The considerations in section 3 will allow us to obtain analytically a two-sided bound on the ratio g 3 /g 2 . 6 We then consider numerically the problem of bounding all coefficients g 2 , f 2 , g 3 , f 3 .

Two-sided analytic bounds
We can derive two-sided bounds analytically using the bracket definitions of g 2 and g 3 , plus the first g-type null constraint, Y 2,0 : We will show that these sum rules lead to a two-sided bound on g 3 /g 2 , Upper bound: Consider the terms entering the brackets in (4.5). First, note that J = J(J + 1), so which implies that Combining these results we find that g 3 is bound by g 2 : showing the upper bound in (4.7).
Lower bound: Now we must include the null constraint (4.6). The problem has a geometric interpretation. We will bound the ratio g 3 /g 2 and we can therefore think of the sum rule for g 2 , (4.4), as a normalization. The right-hand sides of the remaining two sum rules (4.5) and (4.6) then describe a convex hull of points in a two-dimensional space where the allowed values of g 3 /g 2 are given by the intersection of this convex hull with the x-axis. The minimal value (4.7) is found at a point on the boundary of the convex hull, with spectral densities supported only on the odd bracket at spin J = 3 and J = 5. We give more details on this, as well as a graphical visualization, in appendix B.1.
f 2 dependent bound: In fact, we can immediately extend the result just described.
Consider the derivation of the result (4.7) using an arbitrary value of x 1 . The positivity conditions and argument are unchanged, with · · · 0 → · · · x 1 , leading to the more general result r min M 2 (4.12)

Numerical bounds
We will now proceed with a numerical implementation, using optimization with the semidefinite programming solver SDPB [44,45]. We leave a more general discussion for appendix C.2; in the present section we will briefly describe an example with g 2 and g 3 . Let us use the results from (4.4)-(4.6) to define For the numerical implementation we put M 2 = 1 and restore the M 2 dependence in our results using dimensional analysis.
We can now schematically write semidefinite optimization problems of this form  For practical purposes, we limit to a finite set of spins, supplemented by the condition from the limit J → ∞ (see appendix C.2 for more details). Using that g 2 > 0, we can use this algorithm to generate two-sided bounds on g 3 /g 2 . The lower bound follows from minimizing A using the plus sign, to give g 3 /g 2 −A. Conversely, the upper bound follows from choosing the minus sign: giving g 3 /g 2 A (where, of course, A now is different).
With a suitable generalization of the optimization problem (see appendix C.2.1), we are also able to generate bounds involving more than one ratio of variables. In figure 2 we display bounds in the planes spanned by f 2 g 2 , g 3 g 2 and f 3 g 2 , g 3 g 2 . We observe that the bounds on g 3 /g 2 (resp. f 3 /g 2 ) appear to be linear in |f 2 | (resp. f 2 ). A similar observation holds also for the bounds on eight-derivative coefficients considered below. As a consequence, the f 2 dependence in subsequent bounds can be captured by constructing the bound with a few different fixed values of the ratio f 2 /g 2 . In this way, the whole set of bounds involving four-and six-derivative operators can be visualized in a single two-dimensional diagram, which we give in figure 3.
-5 5 10 f 3 /g 2 e a Figure 2. Bounds on the six-derivative terms f 3 /g 2 and g 3 /g 2 as a function of f 2 /g 2 . The dots refer to the values in the partial UV completions discussed in section 4.4: massive axion (a), scalar (φ)and graviton (h), and QED (e), scalar QED (ẽ) and W ± sector (W ).
A special point of interest is the lower bound of g 3 at the special value f 2 = g 2 . We can get a simple analytic value for this point using one null constraint:

Dependence on the number of null constraints
As mentioned before, two-sided bounds can be obtained by using at least one null constraint. On the other hand, we have access to an infinite family of null constraints, and it is interesting to see how our results change when including more of them. First, note that we do not expect any improved results for upper bounds with increased number of null constraints. This can in principle be seen from the explicit form of the null constraints. However, looking ahead, it is clear that it must be the case since the upper bound is saturated by the partial UV completion given by the massive scalar and axion in  Table 1. Lower bounds g 3 /g 2 r min /M 2 using a number N g (N f ) of g-type (f -type) null constraints. compared to a single null constraint is rather modest with respect of our needs, and for the rest of the paper we will use N g + N f = 9 + 9 null constraints.

Bounds on eight-derivative coefficients
In this section we show results that involve coefficients at eight-derivative order. First we will consider the bounds in the (g 4,1 , g 4,2 ) plane. We find that the allowed values of these couplings fall within a triangular shape, whose boundaries we determine analytically. We then consider bounds in the inhomogeneous space involving the coefficient g 3 against different combinations of eight-derivative coefficients.

Bounds involving only eight-derivative coefficients
We start by considering bounds in the plane parametrized by g 4,1 and g 4,2 . The relevant sum rules are (4.17) An explicit analytic consideration in appendix B.2 shows that the allowed values must lie in a region given by a triangle with vertices  . Showing the triangle with the vertices given by (4.18), and, in orange, the bound (4.19) from [4], where allowed region is above the graph.
We can compare with some results from [4] for homogeneous bounds on eight-derivative coefficients. Translated to our parametrization, the bound found in [4] is − 30 7 2g 4,2 g 4,1 + 2g 4,2 6, g 4,1 + 2g 4,2 > 0. Notice that our bound in (4.18) is stronger than (4.19) in the region g 4,1 < 0, but weaker in the region g 4,1 > 0. We then proceed to finding a numerical bound in the (g 4,1 , g 4,2 ) plane. It turns out that the numerical bound is only slightly more constraining than the analytic bound given by the triangle (4.18), and the difference is barely visible. For instance, for g 4,1 /g 2 > 0, even with O(100) null constraints the lower bound on g 4,2 /g 2 is only stronger than the value obtained analytically from (4.18) by at most 10 −5 .

Bounds involving both six-and eight-derivative coefficients
First, we consider the space of couplings (g 4,1 + 2g 4,2 )/g 2 and g 3 /g 2 . We choose this combination of g 4,1 and g 4,2 because it corresponds to a s-t-u symmetric combinations of the amplitudes.
Then, using the algorithm described in appendix C.2.1, we get numerically the allowed region in figure 5. In particular, we consider the extra constraint f 2 /g 2 = 0, ±1 to see how the bounds depend on this ratio. This plot has strong similarities with the figure 8 found in [7] for scalars. In fact, we believe that the amplitudes (specifically, the s − t − u symmetric A[0, 1] amplitude) for the theory saturating the bounds are the same as the amplitudes given in (5.2) of that paper. We know that this is the case for the scalar / axion amplitudes, which sit at the top right, because we have computed them directly. The minimum value g 3 in the red region coincides (after accounting for a factor of 3 due to conventions) with the value of the corresponding region in [7]. Therefore we believe that the theory that lives at the kink in the red region has an A[0, 1] amplitude corresponding to the stu-pole amplitude given in (5.2) in that paper. We do not know what amplitudes lie at the blue kink on the upper left.  Figure 5. Bounds in the plane g 2 M 2 e 1 = g 3 , g 2 M 4 e 2 = g 4,1 + 2g 4,2 normalized to g 2 . The blue, red regions corresponds to f 2 /g 2 = 0, ±1. Intermediate allowed regions can be obtained by linear interpolation.
Using the same strategy, we find a numerical bound in the space (g 3 /g 2 , f 4 /g 2 ) and study how it changes as we vary f 2 . The result is displayed in figure 6.  Figure 6. Allowed region in the plane (g 3 /g 2 , f 4 /g 2 ) for fixed values of f 2 /g 2 = k. The blue, green, red region correspond respectively to k = 0, 0.5, 1. Intermediate values can be obtained by linear interpolation.

Comparison with partial UV completions
In the sections above, we computed rigorous bounds in the EFT coefficients appearing in the low-energy amplitudes (2.6). It may be interesting to compare our bounds with values in some known partial UV completions. We use the word "partial" to emphasize that the theories we have in mind may not themselves be UV-complete theories, but they may be defined in an arbitrarily large energy range M 2 < E 2 < Λ 2 , where Λ is a scale at which further UV completion is necessary. The partial completions we consider consist of integrating out massive fields, either at tree-level or one-loop level. In tables 2 and 3 we give the values of the EFT coefficients following from this procedure. Note that for the case of fields with mass m entering at loop level, we have M 2 = 4m 2 . Any linear combination of the coefficients in tables 2-3 will also be allowed. Table 2. EFT coefficients that result from integrating out various particles of mass M at tree-level.

Tree-level completions
First we consider partial completions which arise from integrating out a single particle at tree-level.

Scalar: Consider the action
The equation of motion for φ reads Plugging that back into the Lagrangian gives Computing the two-to-two photon amplitudes from L andL results in the same set of EFT coefficients, which are recorded in the first line of table 2.
Axion: Starting from the action the low-energy Lagrangian is identical to the scalar with the replacement F F → FF : The values are reported in table 2.
Vector: One can ask if a massive vector provides a similar example of a partial completion.
The answer is no. 7 Consider a massive vector V with field strength G and coupled to the photon through a V AA three-point interaction. One possible choice is However, the amplitudes in this theory do not contain a pole at mass M . This can be understood by noting that the interaction is removable by a field redefinition: which results in the following Lagrangiañ Hence we see that the theory, through the field redefinitions, is equivalent to a noninteracting massive vector plus a single six-derivative interaction. This leads to amplitudes whose behavior is ∼ s 3 at large s, in violation of our assumptions. One might wonder if a different interaction, perhaps with more derivatives, could lead to an interacting partial completion. In fact, this is impossible on general grounds due to the Landau-Yang theorem [46,47]: angular momentum selection rules prevent a massive vector from decaying into two photons.
Graviton: For a massive tensor (see [48] for a review) coupled to photons we have the Lagrangian 28) 7 We thank Callum Jones for helping to clarify this point.
where the kinetic operator is given by In general, we could consider general couplings, where the tensor T takes the form The amplitudes in this theory grow like s 3 at large s. However, the combinationg = − 1 4 g only grows like s 2 , so this is the combination we will consider here. 8 9 For this particular choice, the tensor T then becomes the canonical energy-momentum tensor of F , given by The equations of motion simplify drastically in this case because the energy-momentum tensor is both conserved and traceless. As a result, we can write down the low-energy Lagrangian explicitly: (4.33) The corresponding Wilson coefficients are included in table 2. Let us stress that these values are speculative, since the graviton amplitudes A(s, t) grow like s 2 at fixed u, and therefore marginally violate the Froissart bound.

Loop-level completions
QED provides an important example of a loop-level completion. In this case, the higherderivative coefficients can be extracted from the full one-loop amplitude first computed by Karplus and Neuman [49]. We extracted the coefficients from the expressions given in [24], who report some misprints in earlier literature. The relevant diagram to evaluate is a single 8 Technically we need it to fall off faster than s 2 . The massless graviton also gives s 2 fall-off, but string theory softens the high-energy behavior of the graviton to O(s 2+α u ), with u < 0 in the physical scattering region. So s 2 may obey the bound at infinity, depending on the sign of the corrections, while s 3 has no chance of satisfying the required bound. 9 Only one combination can give the correct high-energy behavior. To see this, consider the field redefi- This removes the interactionghF F but adds a four-and six-derivative self-interaction for F F . Therefore thẽ g coupling can be changed by tuning the four-and six-derivative interactions. However, these interactions give terms in the amplitude that fall off as s 3 , so any deviation from the massive gravity coupling will fall off like s 3 . Table 3. EFT coefficients in loop-level partial UV completions. Here M 2 = 4m 2 i . Recall from (1.4) that a 1 = g 2 +f 2 16 , a 2 = g 2 −f 2 16 .  51,52] [51] [51] electron e running in a loop coupled to the external photons in the standard three-point interaction.
Likewise we may consider scalar QED and vector QED where a massive charged scalar/vector runs in the loop, coupling to the external photons with both three-and four-point interactions -for the relevant diagrams see e.g. [53]. The values for f 2 and g 2 in vector QED agree with the phenomenologically more interesting case of a gauged W boson, "W ± sector", and for this case, as well as scalar QED, the complete one-loop four-point amplitudes were computed in [51].
The bounds we provide in this paper are valid when the theory is infinitely weakly coupled at the scale M , but this is not the case for loop-level completions. The two-loop corrections to QED were computed by Ritus in the case of QED [54] (or [55]) Using α ≈ 1 137 , the corrections in QED are numerically in the order 10 −2 . For our purposes, however, we may consider QED as a partial UV completion in the limit of arbitrarily weak coupling α.

Discussion
In this paper, we have derived a number of two-sided bounds on EFT coefficients appearing in the effective low-energy theory of photons. Our method assumes only that the highenergy amplitudes are analytic on the upper half-plane, bounded at large s, and weakly coupled up to and including the scale s ∼ M 2 . To derive two-sided bounds on the ratios of our coefficients, we have incorporated the use of null constraints that were introduced in [7] for the case of scattering massless scalars.
For one potential application of our bounds, consider a theory with a non-minimally coupled scalar: where the . . . refer to further interactions, which must all have more than four derivatives. If we compute the amplitudes in this theory, we find The higher interactions cannot affect g 2 or g 3 , so we directly conclude that k 1 k 2 0. This can be seen by computing the ratio, M 2 g 3 /g 2 = 1 + 2k 1 k 2 /k 2 1 . Our bounds imply that this ratio must be less than one. Therefore our results, combined with the fact that the massive scalar lives at the boundary of our allowed region, constrain higher-derivative interactions in the massive scalar partial completion.
We have compared our bounds to a number of "partial UV completions" where the photons interact with only a single massive particle. The result was that the massive scalar and massive axion, which coupled to the photon at tree-level, saturated a number of the bounds. The massive graviton and the loop-level completions lie well inside our bounds. The bounds in section 4 display a few other corners that are not populated by known theories. It would be very interesting if physical theories which live at these kinks could be identified -perhaps theories arising from dimensional reduction or braneworld scenarios could be potential candidates. Or perhaps they are not physical, in which case one should ask what assumptions about the S-matrix are required to rule out the amplitudes saturating these bounds.
We also pointed that the putative tree-level completion involving a massive vector is trivial because the interactions may be removed by field redefinitions. This is also the case when considering a graviton with only h µ µ F 2 coupling, as well as the "axi-vector" and "axitensor," which couple to FF . In general, these theories may violate our bounds. This is related to the argument in [57], where it is shown that such a trivial theory -a graviton with only h µ µ F 2 coupling -violates the conjecture of [58] that higher-derivative corrections always increase the Wald entropy of thermodynamically stable black holes. Clearly, such a theory cannot give rise to a unitary S-matrix -the point of [57] was to highlight that the behavior of the S-matrix is a clean, field-redefinition invariant way to diagnose which theories should be taken seriously when making Swampland-like arguments.
There are a number of directions that deserve to be further explored. A crucial one is to understand the effect of loops. The assumption of weak coupling at the scale M means that loop effects of the high-energy theory may be safely ignored, but it limits the scope of the bootstrap method. Because of this assumption, we are only able to use part of the unitarity constraints -that some of ρ i 0 -and not the full constraints, which also bound ρ i from above, because such an upper bound is trivially satisfied at weak coupling, where ρ i 1.
One consequence of the weak coupling assumption is that ρ 5 does not participate in any of the unitarity inequalities we derived in section 3. As a result, we are unable to place any bounds on the coefficients h i . It is interesting that these coefficients are absent in all of the tree-level completions that we studied. We believe a more general unitarity set-up could allow us to bound ρ 5 as well.
Another important direction is to extend the considerations to photons coupled to gravity. The appropriate parametrization of the low-energy amplitudes in that case was outlined in [29]. This has the benefit that it could shed light on the weak gravity conjecture [59,60]. Previous discussions on the use of positivity bounds to prove the WGC include [36,57,61]; so far the so-called t-channel pole appearing from in the amplitude due to graviton exchange has been a major obstacle to a convincing proof. Recently it was shown in [37] that this obstacle can be circumvented in the case of scalar scattering by taking an integral of the subtracted amplitude centered around small t, rather than the t → 0 limit. This method is limited to d > 4, so pursuing this line requires extending our results to d > 4, and/or extending the methods of [37] to d = 4.

A Unitarity
Unitarity requires that all physical states have positive norm. If we have a basis of states labeled by i, j, then unitarity may be expressed by the requirement i|j 0, i.e. that matrix of inner products is positive semi-definite.
We are concerned with photon scattering. The unitarity constraints on amplitudes with spinning external states was nicely reviewed in [38]; our discussion will largely parallel theirs. The basis of states for two-to-two photon scattering is the set of incoming and outgoing two-particle states which transform in irreducible representations of the Poincaré group. Such states are labeled by squared COM energy s, momentum p, spin J, and helicity λ, as well as the helicities of each constituent particle, λ 1 and λ 2 , which equal ±1. Because the particles are indistinguishable, we have | + − = | − + . So we need consider the following basis of states: |a in = |s, p, J, λ, +, + in , |b in = |s, p, J, λ, +, − in , |c in = |s, p, J, λ, −, − in , |a out = |s, p, J, λ, +, + out , |b out = |s, p, J, λ, +, − out , |c out = |s, p, J, λ, −, − out .
The states |a i and |c i exist only for even spins J, while |b i exists for all spins. There is no mixing between different spin eigenstates, so the positive semi-definiteness of the entire matrix of states will imply positive semi-definiteness on the matrix of states at each spin. Furthermore, we have assumed parity invariance so there is no mixing between even and odd parity states. Thus it is convenient to define odd parity eigenstates Note that in last line |3 i is only defined for J ≥ 2. This is because the two-particle state formed from incoming particles with helicity λ 1 and λ 2 is λ = λ 1 − λ 2 . Therefore |3 i has helicity λ = 2. Helicity 2 is only a possible eigenvalue for angular momentum J ≥ 2. Now we can look at the matrix of inner products of all ingoing and outgoing states, and study the corresponding positivity conditions. For the parity odd sector, we find that unitarity implies For parity even states, we have Finally, for even spins greater than 0, the |2 and |3 even parity states can mix, leading to the larger matrix To derive positivity constraints from this, we need where T ij represents the interacting part of the S-matrix. Using this, and applying positive semi-definiteness to each spin individually, we find the following conditions on the partial waves: Positivity of the three-by-three principal minors does not lead to any new positivity conditions. 10 These equations give us positive partial wave expansions for three of our five amplitudes. However, the t-u crossing symmetry which relates A +−−+ and A +−+− leads to another relation: 11 To summarize our positivity properties, we first define the spectral densities ρ i J = Im A i J . These satisfy where we use the weak coupling assumption to simplify the inequalities (A.9)-(A.12). These conditions will allow us to derive sum rules which bound the EFT coefficients defined in (2.6).

B Details on analytic bounds
B.1 Simple analytic bound using one null constraint Here we will derive a two-sided bound on g 3 /g 2 using the problem defined by (4.4)-(4.6). It is instructive to think of these sum rules in terms of a reduced bracket, defined by including the positive quantity 1/(m 4 g 2 ) in the integration measure: Using this definition we can write the sum rules (4.4)-(4.6) as 1 = 1 red.,0 + 1 red.,e + 1 red.,o , (B.2) g 3 g 2 = W 0 g 3 red.,0 + W e g 3 red.,e + W o g 3 red.,o , (B.3) We will now derive the lower bound r min defined by  In figure 7 we plot the vectors W i for J 9 and m 2 M 2 . In the (e 1 , e 2 ) plane, each spin J determines a section of a parabola starting at m 2 = M 2 . The convex hull of all the allowed points for all J is indicated in gray. 12 It is clear that the upper bound comes from the point W 0 (M 2 , 0), giving immediately the bound (4.11). From the figure, we can also see that the lowest possible value of g 3 /g 2 must come from a linear combination of the point Note first that assuming f 2 = g 2 implies that g 3 = 3f 3 , since multiplying (4.12) by g 2 +x 1 f 2 and taking the limit x 1 → −1 shows that 0 g 3 − 1 3 f 3 0. Likewise with f 2 = −g 2 we must have f 3 = −3g 3 .

B.2 Triangular region in the eight-derivative coefficient space
Here we will show how to determine the triangular region in the (g 4,1 , g 4,2 ) plane given by (4.18) in the main text. We will again work with the reduced brackets defined in (B.1), and use the sum rules We then define and consider the following problem: For which A is it possible to find a c 1 such that for all m 2 1 and J in the sums defining the respective brackets? 13 By limiting to spins J J max for some finite J max , the inequality (B.25) can be simplified to a set of inequalities in the (A, c 1 ) plane. The upper and lower value for A resulting from these inequalities does not change by increasing J max as long as J max 5: By finding the intersection with the bound (B.23) we have limited the allowed region in the (g 4,1 , g 4,2 ) plane to the triangle given by (4.18) and shown in figure 4.

C.1 Sum rules
We consider the following system 14 14 For bounds that involve only the constants g k(,i) , for which V + = V − , we will use that Here the sum goes over a range that covers a selection of size N c of the constants c i , which include both the g and f EFT coefficients.
where each entry is a function of m 2 and J, and b = +, −, e, o. This corresponds to writing C g i in a diagonal basis. N f and N g denote the number of f -type and g-type null constraints used.
In table 4 we collect the explicit form of the V b c i for the EFT coefficients considered in this paper. They were derived using the sum rules of section 3.1 in the main text. From eight-derivative order and onwards, they are unique only up to the addition of null constraints.
The null constraints are computed (up to an overall factor) using the replacement rules (3.21)-(3.21) given in the main text. The first two null constraints of each type are included in table 5. Note that we can always choose the form of the null constraints such that

C.2 Optimization problem
Following the previous considerations in section 3, we have a relationship between each low-energy coefficient with a higher energy part of the form where the · b are the positive functional described in (3.8), b = 1, . . . , k and their number depends on whether we are considering f -coefficients or not. We define V b g i (m 2 , J) as the Table 5. Solving the sum rules gives the EFT coefficients in terms of brackets. Below are the terms which appear in each bracket in the solutions. high-energy part of g i corresponding to the k-functional, following the notation in [7]. In the same way we write the null constraints as To get a double-side bound involving the constants g i and g j we follow the same strategy in [7], generalizing it for the case of an arbitrary number of positive functionals · b , b = 1, . . . , k. Let us define  where X represents an arbitrary number of null constraints. With the following algorithm we are able to get a lower and upper bound for g j with respect to g i , or vice versa. First of all, we remark that To obtain an upper bound, we solve the following optimization problem: where c is a set of free coefficients which correspond to the number of null constraints used for defining V . This relation should hold true for each value of spin in the definition of the respective · · · b , J ∈ I b .. If this is the case, it implies that It is easy to see the fundamental importance of positive functionals to relate the high-energy part with the low one, maintaining the bound from the optimization problem. We repeat a similar strategy for the lower bound: In particular, we mainly use g i = g 2 . In this way we can use positivity of g 2 , (4.1) to write a two-sided bound as A g j g 2 B.
(C. 13) In other sections, we use two different versions of this optimization problem. One with 3 positive functionals (3.8), which bounds only g-coefficients and another one with 4 functionals (3.9), which bounds f -and g-coefficients. We will generally usẽ g j := g j g 2 . (C.14)

C.2.1 2d bounds
We can slightly modify the previous optimization problem to get a bound with respect to two low-energy coefficients. In our case, one of them is always g 2 which we take to be positive. Let us define V 1 * (m 2 , J) := V 1 g 2 (m 2 , J), V 1 g i (m 2 , J), V 1 g j (m 2 , J), X 1 (m 2 , J) , V 2 * (m 2 , J) := V 2 g 2 (m 2 , J), V 2 g i (m 2 , J), V 2 g j (m 2 , J), X 2 (m 2 , J) , . . . Using the previous results, we have a lower and upper bound of g i and g l w.r.t g 2 . This implies that if we plot g i and g l , we obtain a rectangular allowed region. We can shrink it solving the following optimization problem where we substitute tog i with values in the finite interval (g i,min ,g i,max ), whereg i,min and g i,max are the values obtained from the previous optimization problem. Mapping many points in the interval, we are able to generate a convex allowed region like in figure 6.

C.2.2 Implementation in SDPB
To solve the optimization problems (C.10) and (C.7) we use SDPB [44,45]. In order to have the right input for the semidefinite optimization problem: • Normalize all V i (m 2 , J) in order to have m 2 only in the numerator. This does not create problems, because it consists in rescaling our system for a positive factor.
• Substitute m 2 → M 2 (1 + x). In this way we have a system of polynomials in x.
We remark that to obtain a stable value we should consider a J max depending on the number of null constraints used. In fact, for computations with ∼ 10 null constraints we need J max ∼ 40, meanwhile, for ∼ 100 null constraints, J max ∼ 120 to get a stable result.