NLO Perturbativity Bounds on Quartic Couplings in Renormalizable Theories with $φ^4$-like Scalar Sectors
Christopher W. Murphy





cmurphy@quark.phy.bnl.gov
Department of Physics, Brookhaven National Laboratory, Upton, N.Y., 11973, U.S.A.


The apparent breakdown of unitarity in low order perturbation theory is often is used to place bounds on the parameters of a theory.
In this work we give an algorithm for approximately computing the next-to-leading order (NLO) perturbativity bounds on the quartic couplings of a renormalizable theory whose scalar sector is ϕ^4-like.
By this we mean theories where either there are no cubic scalar interactions, or the cubic couplings are related to the quartic couplings through spontaneous symmetry breaking.
The quantity that tests where perturbation theory breaks down itself can be written as a perturbative series, and having the NLO terms allows one to test how well the series converges.
We also present a simple example to illustrate the effect of considering these bounds at different orders in perturbation theory.
For example, there is a noticeable difference in the viable parameter when the square of the NLO piece is included versus when it is not.


NLO Perturbativity Bounds on Quartic Couplings in Renormalizable Theories with ϕ^4-like Scalar Sectors
    Christopher W. Murphy
    December 30, 2023
======================================================================================================






§ INTRODUCTION

The unitarity of the S-matrix is frequently used to place theoretical constraints on the parameters of a theory.
On the one hand, if a 2 → 2 scattering amplitude grows with energy, as is typically the case in non-renormalizable theories, then the condition S^† S = 1 will inevitably be violated at some energy scale.
This energy scale then sets an upper limit on where new degrees-of-freedom must appear to unitarize the scattering amplitude. 
While is interesting in its own right, it is not the focus of this work.
On the other hand, in renormalizable theories, there are relations among the parameters of the theory that cancel this growth with energy <cit.>. 
Nevertheless the same procedure can be used to place “perturbativity bounds” on the parameters of a renormalizable theory.
Most famously, a leading order (LO) analysis of this type yielded an upper bound on the mass of the Higgs boson in the Standard Model (SM), m_h  1 TeV <cit.>. 
If some combination of parameters in a renormalizable theory are too large, the amplitude will appear to be non-unitary at some order in perturbation theory. 
Of course these theories are unitary.
The more accurate statement is that perturbation theory is breaking down for this choice of parameters.

This method has subsequently been improved and refined. 
See for instance <cit.> for studies of beyond leading order effects in the SM, including both renormalization group (RG) improvement and higher fixed order contributions.
If some choice of parameters violates the perturbativity bound at tree level it may be that perturbativity is restored at one-loop, or perhaps it may be that those parameters are not viable at any (low) order in perturbation theory.
The preceding discussion calls attention to the following fact.
The quantity that tests where perturbation theory breaks down itself can be written as a perturbative series.
If higher order terms are known, one is then able to test how well the series converges. 
To date not much work has been done along these lines in theories beyond the SM. 
The one-loop corrections necessary to compute perturbativity bounds at next-to-leading order (NLO) in the Two-Higgs Doublet Model (2HDM) with a softly broken 𝐙_2 symmetry were computed in Ref. <cit.>. 
Ref. <cit.> then performed a comprehensive analysis of the viable parameter space in the 2HDM using these NLO perturbativity bounds. 
Prior to this work, an ansatz inspired by SM results was used to to estimate the higher order corrections in the 2HDM <cit.>. 
Not only did the work of Ref.s <cit.> resolve the ambiguity of how to implement perturbativity bounds at NLO in the 2HDM, but they also revealed where the dominate NLO contribution comes from in the SM and the 2HDM.

In this work we first derive a formula for the functional form of the partial wave matrix for high energy 2 → 2 scalar scattering in a general renormalizable theory.
This allows us to construct an algorithm for approximately computing the NLO perturbativity bounds on the scalar quartic couplings of a general renormalizable theory, (approximately) generalizing the result of <cit.>.
Ref. <cit.> showed that this approximation dominates the NLO contribution to the perturbativity bounds in the SM and certain special cases of the 2HDM.
Ref. <cit.> went further and showed that this is generally a good approximation in the 2HDM with a (softly broken) 𝐙_2 symmetry.
Since the approximation is based on the pattern found in <cit.> we expect that this approximation should generally be a good estimate of the full NLO contribution in theories whose scalar sector is ϕ^4-like. 
By this we mean theories where either there are no cubic scalar interactions, or the cubic couplings are related to the quartic couplings through spontaneous symmetry breaking.
 
This approximate NLO result only requires knowledge of the leading order matrix of partial wave amplitudes, and the one-loop scalar contribution to the beta function of each quartic coupling that is to be bounded. 
The advantage of this approximation is its simplicity as both of the required quantities are relatively easy to determine.
To further ease the calculation of these NLO bounds a package implementing this algorithm, 𝙽𝙻𝙾𝚄𝚗𝚒𝚝𝚊𝚛𝚒𝚝𝚢𝙱𝚘𝚞𝚗𝚍𝚜, is available at <https://github.com/christopher-w-murphy/NLOUnitarityBounds> in both Mathematica and Jupyter Notebook formats.

A (likely incomplete) list of other models for which all of the results necessary to implement this algorithm are already known is: the Manohar-Wise Model <cit.>, its extension to include an additional color singlet, SU(2)_L doublet <cit.>, and the Left-Right Symmetric Model <cit.>.

As we will show, the natural interpretation of the couplings appearing in the NLO partial wave amplitudes are RG improved couplings evaluated at a scale much larger than the typical scales of the theory.
As such, these corrections are naturally useful in analyses at high scales. 
Examples of this include high scale flavor alignment <cit.>, effects of custodial symmetry breaking at high energies <cit.>, or simply investigating the validity of a model up to a high energy scale <cit.>.

The rest of this paper starts with a brief review of perturbative unitarity in Sec. <ref>. 
Following that is a derivation of our main results in Sec. <ref>. 
We then present an example to illustrate the effect of considering perturbativity bounds at different orders in perturbation theory in Sec. <ref>.
The example is simple, but it serves to highlight how the viable parameter space in a model can change from order-to-order in perturbation theory.
In particular, there is a noticeable difference in the viable parameter when the square of the NLO piece is included versus when it is not.
The implementation of this model is included in the example notebook associated with the 𝙽𝙻𝙾𝚄𝚗𝚒𝚝𝚊𝚛𝚒𝚝𝚢𝙱𝚘𝚞𝚗𝚍𝚜 package.
Lastly, we discuss our findings in Sec. <ref>.




§ BRIEF REVIEW OF PERTURBATIVE UNITARITY



In this section we give a brief review of perturbative unitarity.  
The S-matrix is unitary, S^† S = 1. 
This condition can be translated into a relation among the various partial wave amplitudes of a given theory, see e.g. <cit.>

    Im(𝐚_j^2 → 2) = (𝐚_j^2 → 2)^†𝐚_j^2 → 2 + ∑_n > 2(𝐚_j^2 → n)^†𝐚_j^2 → n ,

where 𝐚_j^2 → 2 is computed in a basis such that it is diagonal, i.e. it is the eigenvalues of 𝐚_j^2 → 2 that satisfy the relation (<ref>). 
Note also that an integral over the n-body phase space for each term in the sum in the rightmost term in (<ref>) is left understood. 
This is nothing but the equation for an n-sphere of radius 12 centered at Re(𝐚_j^2 → 2) = |𝐚_j^2 → n| = 0, Im(𝐚_j^2 → 2) = 12. 

The imaginary part of 𝐚_j^2 → 2 is fixed by the real part of 𝐚_j^2 → 2 and the 2 → n partial wave amplitudes. 
Focusing on the eigenvalues, a_j^2 → 2, of the matrix 𝐚_j^2 → 2 we have

    2 Im(a_j^2 → 2)_∓ = 1 ∓√(1 - 4 𝒜_j^2) ,   𝒜_j^2 = [Re(a_j^2 → 2)]^2 + ∑_n > 2| a_j^2 → n|^2 ,

with the - (+) solution corresponding to the case when the imaginary part of a_j^2 → 2 is less than (greater than) one-half. 
Assuming a perturbative expansion is viable, the first few orders of Eq. (<ref>) take the form

    Im(a_j^(0))_-    = 0 , 
    Im(a_j^(1))_-    = (a_j^(0))^2 , 
    Im(a_j^(2))_-    = 2 a_j^(0) Re(a_j^(1)) + |a_j^2 → 3, (0)|^2 ,

where the superscript (ℓ) is the perturbative order and we have dropped the superscript 2 → 2. 
For the + solution we instead have Im(a_j^(0))_+ = 1, Im(a_j^(ℓ > 0))_+ = - Im(a_j^(ℓ))_-. 
Assuming a valid perturbative expansion, the + solution corresponds to a scattering amplitude whose tree level imaginary part is Im(ℳ^(0)) = 16 π. 
Since this is not a frequently encountered scenario we will not consider this possibility further in this work, drop the subscript -, and set an upper limit of Im(a_j) ≤12 in the process.

Perturbative unitarity bounds are inequalities derived from Eq. (<ref>). 
From (<ref>) it is clear that the conditions 0 ≤ Im(a_j^2 → 2) ≤12 and 0 ≤𝒜_j^2 ≤14 yield equivalent bounds when the imaginary part of a_j^2 → 2 is kept in its exact form. 
However this is not the case if Im(a_j^2 → 2) is expanded to a finite order in perturbation theory, as in (<ref>). 
For example, to two-loop order, the corresponding bounds are always weaker than those from 0 ≤𝒜_j^2 ≤14. 
Explicitly, to leading order, 0 ≤ Im(a_j^2 → 2) ≈ (a_j^(0))^2 ≤12 whereas 0 ≤𝒜_j^2 ≈ (a_j^(0))^2 ≤14. 
Starting at three-loops in Im(a_j), NNLO for 𝒜_j, the relative strength of the two bounds in no longer fixed.

Once the approximate NLO contributions to the eigenvalues of the partial wave matrix 
 are known, which are derived in Sec. <ref>, perturbative unitarity bounds can be obtained by evaluating one or more of the following

    LO:  (a_0^(0))^2 ≤1/4, 
       NLO:   0 ≤(a_0^(0))^2  + 2 (a_0^(0)) Re(a_0^(1)) ≤1/4 , 
       NLO+:  [(a_0^(0)) + Re(a_0^(1))]^2 ≤1/4.

If perturbation theory is valid it is expected that the bounds obtained from the upper limit (≤14) will be similar in all three cases since the LO eigenvalue is non-trivial and the NLO(+) piece should represent a small correction, but it is important to test, and confirm or deny, if this is actually the case.
The  bound obtained from the lower limit (0 ≤) originates from Ref. <cit.>, Eq. (75) in particular.[This perturbativity bound was subsequently investigated in the 2HDM in Ref. <cit.> where it was called R_1 (not R_1^').] 
In contrast to the ≤14 bound, the 0 ≤ bound only becomes non-trivial at NLO, and can only be violated when perturbation theory breaks down.
While it is expected that the parameter space of a given theory will be more constrained by the NLO perturbativity bounds because of the additional 0 ≤ handle, using this bound goes somewhat against the spirit of the introduction in that there is no similar higher-order term to compare with.
In fact, the NLO+ expression shows how the apparent 0 ≤ violation of unitarity is resolved at higher-orders in perturbation theory. 
(Note that the NLO+ expression contains some NNLO terms, but of course is not the full NNLO expression.) 
To bring things full circle, at NNLO similar apparent 0 ≤ violations of unitarity can occur from the interference between the tree level and two-loop 2 → 2 amplitudes and/or the tree level and one-loop 2 → 3 amplitudes.  
However, there is no reason to expect the higher-order versions of the 0 ≤ bound to be similar to the NLO version, as was the case for the ≤14 bound, since there is little and/or no overlap between the potentially problematic terms.




§ GENERIC PARTIAL WAVE AMPLITUDES FOR HIGH ENERGY SCALAR SCATTERING





 §.§ Elements of the Partial Wave Matrix


Consider a potential of the form

    V = 12 m_α^2 ϕ_αϕ_α + κ_αβγϕ_αϕ_βϕ_γ + λ_αβγδϕ_αϕ_βϕ_γϕ_δ,

where the subscript Greek letters are flavor indices.
The 2 → 2 scattering of high energy (E ≫ |κ_αβγ|, m_α) scalars can schematically by written as

    ℳ_i → f = (Z_ϕ^1/2)^4 𝐕[ϕ^4] ,

where (Z_ϕ^1/2)^4 is the product of the four external wavefunction renormalization factors, and 𝐕[ϕ^4] is the four-point function.
At tree level the unrenormalized four-point function is simply a linear combination of quartic couplings,

    𝐕[ϕ^4]_tree = - c_m λ_m B = - c_m (λ_m + δλ_m) ,

where δλ_m is the counterterm associated with the renormalized coupling λ_m, c_m is the numeric coefficient associated with a given λ_m, and a subscript Roman letter is shorthand for a set of Greek letters, e.g. m = αβγδ.
At high energies the diagrams involving cubic couplings generally do not contribute to Eq. (<ref>).
For s-channel processes this is manifestly true simply because we assume E ≫ |κ_αβγ|.
To see this is also the case for t- and u-channel processes one should retain the full mass dependence of the diagram until after the partial wave amplitude is computed at which point it is safe to take the high energy limit.
An exception to this occurs when all of the particles, internal and external, in a digram are massless.
Such a situation would occur if the theory contains a neutral, CP-even Goldstone boson since it would then be possible to write down vertices with an odd number of this Goldstone boson.
In this case there is a physical divergence in the forward region, analogous to Rutherford scattering, and this method as it is currently implemented is not applicable.
However it is worth pointing that a careful study of the analytic structure of amplitudes involving the t-channel exchange of massless particles showed that the naïve sum rules for processes such as W_L^+ W_L^- scattering are still correct <cit.>, so perhaps there may still be a way to extract perturbativity bounds from such amplitudes. 
In any case, henceforth we will neglect the trilinear couplings κ_αβγ.

A generic one-loop diagram in D = 4 - 2 ε dimensions with four external and two internal scalars, which is the only topology of all scalar, 1PI diagrams that persists in the high energy limit, takes the following form 

    λ_m λ_n/16 π^2(1/ε + 2 - ln(- p^2 - i 0_+/μ^2)) .

As is typically done the scale μ has been introduced to keep the quartic couplings dimensionless. 
The sum of all such diagrams leads to the four-point function at the one-loop level

    𝐕[ϕ^4]    = - c_m λ_m + λ_m λ_n/16 π^2[(σ_mn + τ_mn + υ_mn) (1/ε + 2 + ln(s/μ^2)) . 
       . + i πσ_mn  - τ_mnln(- t/s) - υ_mnln(- u/s) ] ,

where s, t, and u are the usual Mandelstam variables, and the branch cut in the logarithm yields ln(- p^2 - i 0_+) →ln(p^2) - i π for p^2 > 0. 
The one-loop correction is bilinear in the various couplings with the (model dependent) coefficients σ, τ, and υ parameterizing the s-, t-, and u-channel contributions, respectively. 

At one-loop the scalar wavefunction renormalization is finite, which allows the beta function of a quartic coupling to be defined simply as

    β_λ_i = μ∂𝐕[ϕ^4]/∂μ ,

where the particular four-point function entering the definition of the beta function is such that 𝐕[ϕ^4]_tree = - λ_i. 
From this we see that 

    c_m β_λ_m = (σ_mn + τ_mn + υ_mn) λ_m λ_n/8 π^2 ,

which determines the purely scalar contribution to the one-loop beta-function.

After renormalization the scattering amplitude takes the form

    ℳ_i → f   = - c_m λ_m - c_m (δλ_m)_fin. - c_m λ_n (δ Z_mn)_fin. + c_m β_λ_m[1 + ln(√(s)/μ)]  
       - λ_m λ_n/16 π^2[- i πσ_mn + τ_mnln(- t/s) + υ_mnln(- u/s) ] .

The counterterms cancel the divergences arising from the one-loop diagrams, and generically contain finite parts that contribute to the scattering amplitude. 
The diagonal and off-diagonal wavefunction renormalization constants are δ Z_mm and δ Z_mn, respectively, both of which are real (except when using a complex-mass scheme).[In the SM there is a relation between the tree level 2 → 3 partial wave amplitudes and the wavefunction renormalization contribution to the one-loop 2 → 2 partial wave amplitudes <cit.>, which in our notation takes the form

    |a_j^2 → 3, (0)|^2 = δ Z_mm|a_j^2 → 2, (0)|^2   (in the SM).

This causes a partial cancellation in the NLO expression for 𝒜_j that makes our approximation, discussed after (<ref>),  a more accurate representation of 𝒜_j, again at least the in SM. It would be interesting to see if (a generalization of) this relation is true in other theories.] 
The off-diagonal terms generally involve a different linear combination of the tree level amplitudes than the diagonal terms. 
The wavefunction counterterms depend on the particular process under consideration, whereas the δλ_m are process independent.

The full energy dependence of ℳ can be subsumed into a running coupling using standard renormalization group methods

    μ∂λ̅_m(μ)/∂μ = β_λ_m,   λ̅_m(μ_match.) = (λ_m)_phys. ,

with (λ_m)_phys. being the combination of physical parameters that defines λ_m at the scale μ_match..
The scattering amplitude now takes the form

    ℳ_i → f   = - c_m λ̅_m - c_m (δλ_m)_fin. - c_m λ̅_n (δ Z_mn)_fin. + c_m β_λ̅_m
       + λ̅_m λ̅_n/16 π^2[i πσ_mn - τ_mnln(- t/s) - υ_mnln(- u/s) ] .

Retaining only the first term on the right-hand side of Eq. (<ref>) corresponds to the leading-log (LL) approximation of the NLO contribution.

An element of 𝐚_j^2 → 2 is related to the scattering amplitude for the process ℳ_i → f as follows

    (𝐚_j^2 → 2)_i, f = 1/16 π s∫_-s^0  dt  ℳ_i → f(s, t) P_j(1 + 2 ts),

where P_j are the Legendre polynomials. 
In the high energy limit ℳ is independent of s and t at leading order, allowing us to concentrate on the j = 0 case.
Plugging (<ref>) into (<ref>) and simplifying the result using (<ref>) we find

    16 π(𝐚_0^2 → 2)_i,f = - c_m λ̅_m + 3/2 c_m β_λ̅_m + i π - 1/16 π^2σ_mnλ̅_m λ̅_n - c_m (δλ_m)_fin. - c_m λ̅_n (δ Z_mn)_fin.  .

Eq. (<ref>) is the exact expression for the functional form of the partial wave matrix of 2 → 2 scattering amplitudes in a general renormalizable theory in the high energy limit and assuming the scalar quartic couplings are parametrically larger than the gauge and Yukawa couplings.
To the best of our knowledge this result has not been previously been given in the literature.




 §.§ Eigenvalues of the Partial Wave Matrix


In this subsection we give an approximate formula for the NLO corrections to the eigenvalues of the partial wave matrix for high energy scalar scattering. 
One only needs to knows the leading order partial wave matrix and the (scalar one-loop contribution to) the beta functions of the theory under consideration to make use of this approximation.

Recall the well known formula for the NLO perturbations of the eigenvalues of an eigensystem that is known completely at LO

    a_0^(1) = x⃗^⊤_(0)·𝐚_0^(1)·x⃗_(0) .

It says that the NLO eigenvalues depend only on the NLO correction to the matrix and the LO eigenvectors, which are determined from the LO matrix. 
Here the exact and leading order eigensystems respectively are

    𝐚_0 ·x⃗   = a_0  x⃗ , 
    𝐚_0^(0)·x⃗_(0)   = a_0^(0) x⃗_(0) ,

and each object appearing in the first line of (<ref>) is assumed to have an expansion

    𝐚_0    = 𝐚_0^(0) + 𝐚_0^(1) + … , 
    x⃗   = x⃗_(0) + x⃗_(1) + … , 
    
    a_0    = a_0^(0) + a_0^(1) + … .


The second term on the right-hand side of Eq. (<ref>) (the β-function contribution) to an element of 𝐚_0 is universal. 
Therefore the β-function contributions to the eigenvalues of 𝐚_0 can simply be determined using Eq (<ref>). 
All one has to do to obtain 𝐚_0,β^(1) from 𝐚_0^(0) is replace each quartic coupling λ that appears in 𝐚_0^(0) with - 32β_λ. 
Then using the LO eigenvalues x⃗_(0), which are also determined from 𝐚_0^(0), we find the β contribution to a_0^(1)

    a_0, β^(1)   =  x⃗_(0)^⊤·𝐚_0, β^(1)·x⃗_(0) , 
    𝐚_0,β^(1)   = - 3/2. 𝐚_0^(0)|_λ_m →β_λ_m .


Unlike the β-function contribution, the third term on the right-hand side of Eq. (<ref>) (the σ contribution) is not universal.
Without doing an explicit calculation, the σ term contribution to the partial wave matrix is not known.
However we are interested in knowing the eigenvalues of the matrix rather than the entire matrix itself.
Knowing that the imaginary parts of the NLO elements of 𝐚_0 come exclusively from the σ contribution, we can use the second line of (<ref>), i.e. the fact that the theory is unitary, to determine the σ contribution to the NLO eigenvalues. 
Note that the numerical factor in front of σ_mn in Eq. (<ref>) is complex.
This is convenient as it allows us to use the unitarity of the theory to also determine the real part of the sigma contribution the NLO eigenvalues.
This yields the final result for the sigma contribution

    a_0, σ^(1) =  (i - 1/π) (a_0^(0))^2 .
  

To proceed further without doing an explicit calculation we must drop the finite pieces of the counterterms.
Thus combining the β contribution and the σ contribution we arrive at our result for the approximation NLO contribution to the eigenvalues of the partial wave matrix

    a_0^(1) = a_0,β^(1) + a_0,σ^(1) .





§ A SIMPLE EXAMPLE



In this section we present a simple example to sketch the effect of using perturbativity bounds at different orders in perturbation theory.
Consider as a toy model the 2HDM with a U(2) symmetry, instead of 𝐙_2, to prevent tree level flavor changing neutral currents <cit.>.
The scalar sector of this model has only two unique quartic couplings, and the U(2) symmetry preserves this relation along the RG-flow.
The potential is

    V = λ_1/2[(H_1^† H_1) + (H_2^† H_2) - v^2/2]^2  + (λ_1 - λ_3) [(H_1^† H_2) (H_2^† H_1) - (H_1^† H_1) (H_2^† H_2)],

where H_1,2 are the two Higgs doublets.[One could simplify the potential in this case by defining λ_4 ≡λ_1 - λ_3, but we stick with λ_1 and λ_3 as they are what is used in <cit.>.]
We will take λ_1,3 to be free parameters for the purposes of this demonstration, and neglect the electroweak vev, v ≈ 246 GeV, as it is not important at high energies.

Figure <ref> shows the bounds on λ_1 and λ_3 that result when the various perturbativity conditions, (<ref>), are applied. 
The parameter space shaded blue, orange, green, red, and purple is ruled out by perturbativity bounds from eigenvalues whose LO terms are proportional to 4 λ_1 + λ_3, - λ_1 + 2 λ_3, 2 λ_1 - λ_3, λ_1, and λ_3, respectively. 
There is a noticeable  difference in the viable parameter when the square of the NLO piece is included versus when it is not, as shown in panels (b) and (c) of Fig. <ref>. 
This is reminiscent of a result in Ref. <cit.>, where the bounds obtained on the coefficients of some dimension-6 operators in the SMEFT depend noticeably on whether or not the square of the dimension-6 amplitude is included in the calculation of the cross section under consideration. 
This parameter space could be further constrained theoretically by required the potential to be bounded from below, etc.



Figure <ref> shows in blue, green, and red the viable parameter space for two individual eigenvalues based on the LO, NLO, and NLO+ perturbativity conditions, respectively. 
The contours in panels (a) and (b) of Fig. <ref> are colored blue and orange, respectively, to match the color coding of the eigenvalues in Fig. <ref>.
The parameter space viable at NLO (green) in Fig. <ref> is determined using both the upper and lower limits on a_0, which, in contrast, are illustrated separately in Fig. <ref> in panels (b) and (d), respectively.


The different perturbativity criteria yield different viable parameter spaces.
Perhaps the simplest way to test the convergence of a_0 as a perturbative series is to see which choices of parameters are considered viable by multiple perturbativity criteria.




§ DISCUSSION



In the section we discuss computing the leading order partial wave matrix and one-loop scalar beta function, and the validity of our approximation, before summarizing our algorithm for finding the NLO perturbativity bounds in a given theory.

One of the advantages of this approach is its simplicity in that it only relies on knowledge of the leading order partial wave matrix and one-loop scalar beta function.
In particular, the renormalization group equations in a general quantum field theory to two-loop order have been known for some time, and software exists to derive them in a specific model <cit.>. 
In the extreme energy limit s ≫ M_i^2, which is sufficient to consider if one is only interested in bounding dimensionless couplings, the method of Ginzburg and Ivanov can be used to simply the computation of 𝐚_0^(0) <cit.>. 
Though they originally considered the 2HDM, their argument can be used for any renormalizable theory where the matter fields have definite representations under the gauge group of theory. 
In particular at high energies masses and mixings are not important, and 𝐚_0^(0) will be block diagonal with blocks of definite representations under the gauge and global symmetries of the theory <cit.>. 
In gauge theories, scattering amplitudes involving longitudinally polarized vector bosons should be included when determining the bounds on the quartic couplings of the theory. 
Their inclusion can be greatly simplified through the use of the Goldstone Boson Equivalence Theorem, see e.g. <cit.>. 
Lastly, we note that in gauge theories (with spontaneous symmetry breaking) there is the additional complication of preserving gauge invariance. 
In a mixed MS/on-shell scheme, gauge independence can be spoiled unless tadpole diagrams are properly taken into account. 
This can be done by generalizing the SM results of Fleischer and Jegerlehner <cit.> to the theory under consideration. 
In fact this has been done in the 2HDM <cit.> and the SM Effective Field Theory (SMEFT) <cit.>.

In considering how well our approximation does at capturing the exact NLO results, theories can be placed into one of three categories, broadly speaking.
The first class is theories are those without spontaneous symmetry breaking. 
Here there exists a renormalization scheme such that the approximation is actually exact.
This is simply because all of the counterterms for the quartic couplings are independent of the mass counterterms, and the finite parts of the quartic couplings can be chosen to cancel any potential wavefunction counterterm contribution.
The second class of theories is those with spontaneous symmetry breaking and a ϕ^4-like scalar sector.
This class includes the SM and the 2HDM, theories for which this approximation is known to be good one <cit.>.
In multi-ϕ^4 theory with spontaneous symmetry breaking the wavefunction renormalization is finite at one-loop, but more importantly it is parametrically identical to the 1PI one-loop contribution.
Thus, in ϕ^4-like theories is that there is no qualitative difference between the different NLO contributions, and there is a qualitative difference between the LO and NLO contributions
This suggests that at the very least including some NLO contributions should give a good qualitative estimate of the full NLO contribution, which may be sufficient if one wishes to test the convergence of a_0 in perturbation theory through NLO.
Finally, the third class of theories are those with spontaneous symmetry breaking and scalar cubic interactions whose cubic couplings are not related to its quartic couplings.
In contrast with the first two classes of theories, it is not necessarily the case that this approximation will give a good description of the exact NLO result.
The reason being that the wavefunction renormalization contribution is parametrically different from both the LO and the 1PI NLO contributions to the partial wave amplitude.
On the other hand, if the cubic coupling of interest is much smaller than the internal masses in the wavefunction renormalization diagrams, the approximation may still be a good one as the theory is approaching the ϕ^4 limit in this scenario.
Examples of theories of this type include extending the SM with a real scalar singlet without a 𝐙_2 symmetry, or with scalars that are color singlets and SU(2)_L triplets. This includes the original Georgi-Machacek model, but not its generalizations <cit.>. Additionally, the 2HDM plus a pseudoscalar singlet falls into this class, see for example <cit.> and the references therein.
Another thing to keep in mind in our approach to finding the NLO perturbativity bounds is the quartic couplings entering into the partial-wave amplitudes are necessarily running couplings evaluated at an energy scale much larger than the other scales in the problem. 
A LO analysis is simpler in the sense that expressions could simply involve ordinary quartic couplings. 
One could always RG-improve the LO bounds by replacing the ordinary couplings with running couplings with no penalty.
This is the leading-log approximation. 
However, to do the opposite, replace running couplings in the NLO expressions with ordinary couplings, would be a further approximation.

To summarize, in this work we first derive a formula for the functional form of the partial wave matrix for high energy 2 → 2 scalar scattering in a general renormalizable theory.
This allows us to construct an algorithm for approximately computing the NLO perturbativity bounds on the scalar quartic couplings of a general renormalizable theory.
We expect the approximation to be a good estimate of the full NLO contribution in theories whose scalar sector is ϕ^4-like. 
By this we mean theories where either there are no cubic scalar interactions, or the cubic couplings are related to the quartic couplings through spontaneous symmetry breaking.
The approximate NLO result only requires knowledge of the leading order matrix of partial wave amplitudes, and the one-loop scalar contribution to the beta function of each quartic coupling that is to be bounded, both of which are quantities that are relatively easy to determine.
The algorithm for finding the eigenvalues of the partial wave matrix at approximate next-to-leading order is as follows:


  * Given the leading order partial wave matrix, 𝐚_0^(0), find its eigenvalues, a_0^(0), and eigenvectors, x⃗_(0). 
This may need to be done numerically.

  * The β-function contribution to the NLO eigenvalues, a_0,β^(1), is given by Eq. (<ref>) with

    𝐚_0,β^(1) = - 3/2. 𝐚_0^(0)|_λ_m →β_λ_m .


  * Using the unitarity of the theory, the σ-contribution to the NLO eigenvalues is

    a_0,σ^(1) = (i - 1/π) (a_0^(0))^2 .


  * The NLO contribution is given by the sum of the two pieces, a_0^(1) = a_0,β^(1) + a_0,σ^(1). 

In addition, Mathematica and Jupyter Notebook packages implementing this algorithm, 𝙽𝙻𝙾𝚄𝚗𝚒𝚝𝚊𝚛𝚒𝚝𝚢𝙱𝚘𝚞𝚗𝚍𝚜, are available at <https://github.com/christopher-w-murphy/NLOUnitarityBounds>.
The natural interpretation of the couplings appearing in the NLO partial wave amplitudes are RG improved couplings evaluated at a scale much larger than the typical scales of the theory.
Therefore these corrections are naturally useful in analyses at high scales. 
Once a_0^(1) is known perturbativity bounds on quartic couplings can be obtained at either leading order or next-to-leading order in theories with ϕ^4-like scalar sectors, enabling a test of the convergence of the a_0 as a perturbative series among other studies.
This is highlighted in the example we present.
Specifically, there is a noticeable difference in the viable parameter when the square of the NLO piece is included versus when it is not.



We thank D. Chowdhury, H. Davoudiasl, S. Dawson, O. Eberhardt, B. Grinstein, and P. Uttayarat for useful discussions. This work was supported by the US DOE under grant contract DE-SC0012704.



10

Gunion:1990kf
J. F. Gunion, H. E. Haber, and J. Wudka, “Sum rules for Higgs bosons,”
http://dx.doi.org/10.1103/PhysRevD.43.904Phys. Rev. D43 (1991) 904–912.


Dicus:1992vj
D. A. Dicus and V. S. Mathur, “Upper bounds on the values of masses in
  unified gauge theories,”
http://dx.doi.org/10.1103/PhysRevD.7.3111Phys. Rev. D7 (1973) 3111–3114.


Lee:1977yc
B. W. Lee, C. Quigg, and H. B. Thacker, “The Strength of Weak Interactions at
  Very High-Energies and the Higgs Boson Mass,”
http://dx.doi.org/10.1103/PhysRevLett.38.883Phys. Rev. Lett.
  38 (1977) 883–885.


Lee:1977eg
B. W. Lee, C. Quigg, and H. B. Thacker, “Weak Interactions at Very
  High-Energies: The Role of the Higgs Boson Mass,”
http://dx.doi.org/10.1103/PhysRevD.16.1519Phys. Rev. D16 (1977) 1519.


Marciano:1989ns
W. J. Marciano, G. Valencia, and S. Willenbrock, “Renormalization Group
  Improved Unitarity Bounds on the Higgs Boson and Top Quark Masses,”
http://dx.doi.org/10.1103/PhysRevD.40.1725Phys. Rev. D40 (1989) 1725.


Dawson:1988va
S. Dawson and S. Willenbrock, “Unitarity Constraints on Heavy Higgs
  Bosons,”
http://dx.doi.org/10.1103/PhysRevLett.62.1232Phys. Rev. Lett.
  62 (1989) 1232.


Durand:1992wb
L. Durand, J. M. Johnson, and J. L. Lopez, “Perturbative unitarity and
  high-energy W(L)+-, Z(L), H scattering. One loop corrections and the Higgs
  boson coupling,”
http://dx.doi.org/10.1103/PhysRevD.45.3112Phys. Rev. D45 (1992) 3112–3127.


Durand:1993vn
L. Durand, P. N. Maher, and K. Riesselmann, “Two loop unitarity constraints
  on the Higgs boson coupling,”
  http://dx.doi.org/10.1103/PhysRevD.48.1084Phys. Rev. D48 (1993) 1084–1096,
http://arxiv.org/abs/hep-ph/9303234arXiv:hep-ph/9303234
  [hep-ph].


Grinstein:2015rtl
B. Grinstein, C. W. Murphy, and P. Uttayarat, “One-loop corrections to the
  perturbative unitarity bounds in the CP-conserving two-Higgs doublet model
  with a softly broken ℤ_2 symmetry,”
  http://dx.doi.org/10.1007/JHEP06(2016)070JHEP 06
  (2016) 070,
http://arxiv.org/abs/1512.04567arXiv:1512.04567 [hep-ph].


Cacchio:2016qyh
V. Cacchio, D. Chowdhury, O. Eberhardt, and C. W. Murphy, “Next-to-leading
  order unitarity fits in Two-Higgs-Doublet models with soft ℤ_2
  breaking,” http://dx.doi.org/10.1007/JHEP11(2016)026JHEP
  11 (2016) 026,
http://arxiv.org/abs/1609.01290arXiv:1609.01290 [hep-ph].


Baglio:2014nea
J. Baglio, O. Eberhardt, U. Nierste, and M. Wiebusch, “Benchmarks for Higgs
  Pair Production and Heavy Higgs boson Searches in the Two-Higgs-Doublet Model
  of Type II,” http://dx.doi.org/10.1103/PhysRevD.90.015008Phys.
  Rev. D90 no. 1, (2014) 015008,
http://arxiv.org/abs/1403.1264arXiv:1403.1264 [hep-ph].


Chowdhury:2015yja
D. Chowdhury and O. Eberhardt, “Global fits of the two-loop renormalized
  Two-Higgs-Doublet model with soft Z_2 breaking,”
  http://dx.doi.org/10.1007/JHEP11(2015)052JHEP 11
  (2015) 052,
http://arxiv.org/abs/1503.08216arXiv:1503.08216 [hep-ph].


He:2013tla
X.-G. He, H. Phoon, Y. Tang, and G. Valencia, “Unitarity and vacuum stability
  constraints on the couplings of color octet scalars,”
  http://dx.doi.org/10.1007/JHEP05(2013)026JHEP 05
  (2013) 026,
http://arxiv.org/abs/1303.4848arXiv:1303.4848 [hep-ph].


Cheng:2016tlc
L. Cheng and G. Valencia, “Two Higgs doublet models augmented by a scalar
  colour octet,” http://dx.doi.org/10.1007/JHEP09(2016)079JHEP
  09 (2016) 079,
http://arxiv.org/abs/1606.01298arXiv:1606.01298 [hep-ph].


Cheng:2017tbn
L. Cheng and G. Valencia, “Validity of two Higgs doublet models with a scalar
  color octet up to a high energy scale,”
http://arxiv.org/abs/1703.03445arXiv:1703.03445 [hep-ph].


Chakrabortty:2016wkl
J. Chakrabortty, J. Gluza, T. Jelinski, and T. Srivastava, “Theoretical
  constraints on masses of heavy particles in Left-Right Symmetric Models,”
  http://dx.doi.org/10.1016/j.physletb.2016.05.092Phys. Lett.
  B759 (2016) 361–368,
http://arxiv.org/abs/1604.06987arXiv:1604.06987 [hep-ph].


Gori:2017qwg
S. Gori, H. E. Haber, and E. Santos, “High scale flavor alignment in
  two-Higgs doublet models and its phenomenology,”
http://arxiv.org/abs/1703.05873arXiv:1703.05873 [hep-ph].


Blasi:2017xmc
S. Blasi, S. De Curtis, and K. Yagyu, “Effects of custodial symmetry breaking
  in the Georgi-Machacek model at high energies,”
http://arxiv.org/abs/1704.08512arXiv:1704.08512 [hep-ph].


Grinstein:2013npa
B. Grinstein and P. Uttayarat, “Carving Out Parameter Space in Type-II Two
  Higgs Doublets Model,” http://dx.doi.org/10.1007/JHEP09(2013)110,
  10.1007/JHEP06(2013)094JHEP 06 (2013) 094,
  http://arxiv.org/abs/1304.0028arXiv:1304.0028 [hep-ph].
[Erratum: JHEP09,110(2013)].


Chakrabarty:2014aya
N. Chakrabarty, U. K. Dey, and B. Mukhopadhyaya, “High-scale validity of a
  two-Higgs doublet scenario: a study including LHC data,”
  http://dx.doi.org/10.1007/JHEP12(2014)166JHEP 12
  (2014) 166,
http://arxiv.org/abs/1407.2145arXiv:1407.2145 [hep-ph].


Ferreira:2015rha
P. Ferreira, H. E. Haber, and E. Santos, “Preserving the validity of the
  Two-Higgs Doublet Model up to the Planck scale,”
  http://dx.doi.org/10.1103/PhysRevD.92.033003,
  10.1103/PhysRevD.94.059903Phys. Rev. D92 (2015) 033003,
  http://arxiv.org/abs/1505.04001arXiv:1505.04001
  [hep-ph].
[Erratum: Phys. Rev.D94,no.5,059903(2016)].


Itzykson:1980rh
C. Itzykson and J. B. Zuber, Quantum Field Theory.
International Series In Pure and Applied Physics. McGraw-Hill, New
  York, 1980.
<http://dx.doi.org/10.1063/1.2916419>.


Kilian:2014zja
W. Kilian, T. Ohl, J. Reuter, and M. Sekulla, “High-Energy Vector Boson
  Scattering after the Higgs Discovery,”
  http://dx.doi.org/10.1103/PhysRevD.91.096007Phys. Rev.
  D91 (2015) 096007,
http://arxiv.org/abs/1408.6207arXiv:1408.6207 [hep-ph].


Bellazzini:2014waa
B. Bellazzini, L. Martucci, and R. Torre, “Symmetries, Sum Rules and
  Constraints on Effective Field Theories,”
  http://dx.doi.org/10.1007/JHEP09(2014)100JHEP 09
  (2014) 100,
http://arxiv.org/abs/1405.2960arXiv:1405.2960 [hep-th].


Ivanov:2006yq
I. P. Ivanov, “Minkowski space structure of the Higgs potential in 2HDM,”
  http://dx.doi.org/10.1103/PhysRevD.76.039902,
  10.1103/PhysRevD.75.035001Phys. Rev. D75 (2007) 035001,
  http://arxiv.org/abs/hep-ph/0609018arXiv:hep-ph/0609018
  [hep-ph].
[Erratum: Phys. Rev.D76,039902(2007)].


Ferreira:2009wh
P. M. Ferreira, H. E. Haber, and J. P. Silva, “Generalized CP symmetries and
  special regions of parameter space in the two-Higgs-doublet model,”
  http://dx.doi.org/10.1103/PhysRevD.79.116004Phys. Rev.
  D79 (2009) 116004,
http://arxiv.org/abs/0902.1537arXiv:0902.1537 [hep-ph].


Branco:2011iw
G. C. Branco, P. M. Ferreira, L. Lavoura, M. N. Rebelo, M. Sher, and J. P.
  Silva, “Theory and phenomenology of two-Higgs-doublet models,”
  http://dx.doi.org/10.1016/j.physrep.2012.02.002Phys. Rept.
  516 (2012) 1–102,
http://arxiv.org/abs/1106.0034arXiv:1106.0034 [hep-ph].


Falkowski:2016cxu
A. Falkowski, M. Gonzalez-Alonso, A. Greljo, D. Marzocca, and M. Son,
  “Anomalous Triple Gauge Couplings in the Effective Field Theory Approach at
  the LHC,” http://dx.doi.org/10.1007/JHEP02(2017)115JHEP
  02 (2017) 115,
http://arxiv.org/abs/1609.06312arXiv:1609.06312 [hep-ph].


Lyonnet:2013dna
F. Lyonnet, I. Schienbein, F. Staub, and A. Wingerter, “PyR@TE:
  Renormalization Group Equations for General Gauge Theories,”
  http://dx.doi.org/10.1016/j.cpc.2013.12.002Comput. Phys.
  Commun. 185 (2014) 1130–1152,
http://arxiv.org/abs/1309.7030arXiv:1309.7030 [hep-ph].


Ginzburg:2005dt
I. F. Ginzburg and I. P. Ivanov, “Tree-level unitarity constraints in the
  most general 2HDM,”
  http://dx.doi.org/10.1103/PhysRevD.72.115010Phys. Rev.
  D72 (2005) 115010,
http://arxiv.org/abs/hep-ph/0508020arXiv:hep-ph/0508020
  [hep-ph].


Bagger:1989fc
J. Bagger and C. Schmidt, “Equivalence Theorem Redux,”
http://dx.doi.org/10.1103/PhysRevD.41.264Phys. Rev. D41 (1990) 264.


Fleischer:1980ub
J. Fleischer and F. Jegerlehner, “Radiative Corrections to Higgs Decays in
  the Extended Weinberg-Salam Model,”
http://dx.doi.org/10.1103/PhysRevD.23.2001Phys. Rev. D23 (1981) 2001–2026.


Krause:2016oke
M. Krause, R. Lorenz, M. Muhlleitner, R. Santos, and H. Ziesche,
  “Gauge-independent Renormalization of the 2-Higgs-Doublet Model,”
  http://dx.doi.org/10.1007/JHEP09(2016)143JHEP 09
  (2016) 143,
http://arxiv.org/abs/1605.04853arXiv:1605.04853 [hep-ph].


Denner:2016etu
A. Denner, L. Jenniches, J.-N. Lang, and C. Sturm, “Gauge-independent
  MS renormalization in the 2HDM,”
  http://dx.doi.org/10.1007/JHEP09(2016)115JHEP 09
  (2016) 115,
http://arxiv.org/abs/1607.07352arXiv:1607.07352 [hep-ph].


Krause:2016xku
M. Krause, M. Muhlleitner, R. Santos, and H. Ziesche, “Higgs-to-Higgs boson
  decays in a 2HDM at next-to-leading order,”
  http://dx.doi.org/10.1103/PhysRevD.95.075019Phys. Rev.
  D95 no. 7, (2017) 075019,
http://arxiv.org/abs/1609.04185arXiv:1609.04185 [hep-ph].


Hartmann:2016pil
C. Hartmann, W. Shepherd, and M. Trott, “The Z decay width in the SMEFT:
  y_t and λ corrections at one loop,”
  http://dx.doi.org/10.1007/JHEP03(2017)060JHEP 03
  (2017) 060,
http://arxiv.org/abs/1611.09879arXiv:1611.09879 [hep-ph].


Logan:2015xpa
H. E. Logan and V. Rentala, “All the generalized Georgi-Machacek models,”
  http://dx.doi.org/10.1103/PhysRevD.92.075011Phys. Rev.
  D92 no. 7, (2015) 075011,
http://arxiv.org/abs/1502.01275arXiv:1502.01275 [hep-ph].


Goncalves:2016iyg
D. Goncalves, P. A. N. Machado, and J. M. No, “Simplified Models for Dark
  Matter Face their Consistent Completions,”
  http://dx.doi.org/10.1103/PhysRevD.95.055027Phys. Rev.
  D95 no. 5, (2017) 055027,
http://arxiv.org/abs/1611.04593arXiv:1611.04593 [hep-ph].






