A global analysis of axion-like particle interactions using SMEFT fits

In the presence of an axion or axion-like particle (ALP) that couples to the Standard Model via dimension-five interactions, dimension-six SMEFT interactions are generated via renormalization-group evolution. As many of these SMEFT contributions are experimentally tightly constrained, this"ALP-SMEFT interference"can be used to derive indirect bounds on the ALP couplings to the Standard Model particles. We present a global analysis of the Wilson coefficients of the ALP effective Lagrangian based on Higgs, top, and low-energy data. The obtained bounds are model independent and are competitive or even stronger than direct bounds in the GeV to TeV ALP-mass range.


Introduction
The lack of new-particle discoveries at the Large Hadron Collider (LHC) other than the Higgs boson restricts the space of possibilities for physics beyond the Standard Model (SM). To have evaded detection so far, new particles must either be too heavy to be produced in current highenergy experiments, or interact very weakly with the SM particles. Pseudo Nambu-Goldstone bosons that appear from the spontaneous breaking of a global symmetry, commonly referred to as axions and axion-like particles (ALPs), are among the best motivated light particles beyond the SM. They could play a crucial role in the solution to the strong CP problem [1][2][3] and are potential dark matter candidates (see [4] for a review). While most explicit models addressing the strong CP problem predict a strict relation between the axion mass and its decay constant, it is also possible to obtain solutions to the strong CP problem with heavier axions [5][6][7][8][9][10][11][12][13][14][15][16][17]. Furthermore, composite-Higgs models [18,19] and even supersymmetric extensions of the SM [20][21][22] naturally give rise to ALPs. These models provide a good motivation to search for light pseudoscalar particles in a parameter space that extends beyond the traditional axion searches. Indeed, one of the main interests in searching for ALPs is that they can be a forerunner of a high-scale new physics sector, which otherwise could be hard to access experimentally.
Direct ALP searches span from cosmological [23,24] and astrophysical [25,26] observations to collider [27][28][29][30][31][32] and flavor [33][34][35][36][37][38] experiments. These observables impose powerful constraints on the ALP parameter space, although often under very specific assumptions on its interactions. For example, direct searches for ALPs at particle colliders must make an assumption about the process in which the ALP is produced (e.g. in the decay of a Higgs boson, and Z boson, or in weak decays of kaons or B mesons), about its lifetime, and the way in which it decays, which may involve a decay outside of the detector. In the derivation of astrophysical bounds on, e.g., the ALP-photon coupling from supernovae or ALPs produced in the sun, it makes a crucial difference whether or not a vanishing ALP-electron coupling is assumed. In the vast majority of the existing analyses, it was assumed that only a single ALP couplings is non-zero at the scale of Peccei-Quinn symmetry breaking -an assumption which is not valid in even the classic KSVZ [39,40] and DFSZ [41,42] models for the QCD axion.
As recently shown in [43], the presence of ALP couplings to SM particles yields a non-trivial renormalization-group (RG) flow into the Wilson coefficients of the dimension-six Standard Model Effective Field Theory (SMEFT) operators. This opens up the interesting possibility of indirectly testing ALP interactions by utilizing existing SMEFT studies. In this paper, we exploit the ALP -SMEFT interference to deduce constraints from a global fit including low-energy, Higgs and top data, as implemented in current global SMEFT analyses [44][45][46][47][48][49][50][51][52][53][54][55][56][57][58][59][60][61][62], thus providing a complementary use of these observables to constrain light new physics. As we show, our analysis probes previously unexplored parts of the ALP parameter space, especially for ALP masses above 10 GeV. Furthermore, one of the significant advantages of our indirect approach is that the derived bounds are mostly independent of the ALP mass and require no assumptions on other ALP couplings or the ALP lifetime. This feature overcomes one of the main limitations of direct ALP searches, which often produce highly model-dependent constraints. Therefore, the indirect ALP bounds presented in this paper offer an important complementary approach with respect to direct constraints even when the latter appear to be stronger under certain conditions. This paper is organized as follows: Section 2 introduces the ALP -SMEFT Lagrangian and discusses its connection to the SMEFT. In Section 3, we describe the global analysis setup, present the fit results and compare them to bounds from direct searches. Section 4 is dedicated to a reinterpretation of our global analysis in terms of concrete ultraviolet (UV) completions. We conclude in Section 5.

ALP couplings to the SM
We consider an extension of the SM that includes a pseudoscalar, gauge-singlet ALP as an additional light state. Its couplings to SM fields are, at the classical level, protected by an approximate shift symmetry a → a+c, broken softly by the ALP mass term m a . We relegate the discussion of possible UV completions of this model to Section 4.

The ALP Lagrangian
It has been pointed out in [43] that a consistent effective field theory (EFT) for an extension of the SM featuring an ALP must necessarily include higher-dimensional operators built out of SM fields only. Above the electroweak scale, the most general effective Lagrangian describing the interactions of an ALP with SM particles thus reads with L SMEFT denoting the SMEFT Lagrangian. The SMEFT Lagrangian up to dimension-six order can be found in [63]. In what follows, we adopt the same conventions as in this reference, except for the labels for the Higgs field and its (tachyonic) mass, for which we use H and m H , respectively. Note, in particular, that we work with dimensionless Wilson coefficients and factor out inverse powers of the new-physics scale Λ when writing down higher-dimensional operators in the effective Lagrangian. The Lagrangian L SM+ALP describes the ALP interactions with SM particles. Apart from the soft symmetry breaking by the ALP mass m a , we consider only classically shift-invariant interactions with the SM fields. Specifically, the ALP Lagrangian up to dimension-six operators is then given by 1 where c F (F = q, u, d, l, e) are 3 × 3 hermitian matrices and c ii (i = G, W, B, H) are real couplings.
The real parameter f is referred to as the ALP decay constant, which is related to the relevant new-physics scale by The Lagrangian above can be cast into an alternative but equivalent form, in which the ALP couplings to fermions are of non-derivative type. This is achieved by means of the field redefinitions (4) Note that due to the axial anomaly, additional contributions enter the ALP couplings to gauge bosons, such that c V V and C V V are related by [64] with N c = 3 and the hypercharges Y u = 2/3, The new ALP-fermion coupling matrices are related to the original Lagrangian couplings bỹ with F R = u, d, e and F L = q (l) for quark (lepton) couplings. Assuming a flavor-universal structure for the c F couplings, i.e. c F = c F 1 3 , these interactions reduce tõ where C u,d ≡ c u,d − c q and C e ≡ c e − c l . In the flavor-universal scenario, dimension-five ALP interactions with SM particles are thus fully described by six free parameters: C GG , C W W , C BB , C u , C d and C e . 2

ALP -SMEFT interference
The dimension-five ALP couplings to SM particles generate divergent amplitudes that require dimension-six SMEFT counterterms. As a result, the RG equations for the SMEFT operators get modified by additional inhomogeneous source terms, such that the scale dependence in the presence of the ALP is given by Here γ SMEFT denotes the anomalous-dimension matrix of the dimension-six SMEFT operators calculated in [65][66][67][68], and the source terms γ SMEFT−ALP have been calculated in [43]. 3 The presence of ALP interactions with the SM particles thus generates a non-trivial RG flow into the SMEFT Wilson coefficients. The dimension-six inhomogeneous source terms are independent of the ALP mass. Thus, as long as the ALP mass lies below the scale of the observables used in the fit, we will obtain mass-independent indirect bounds on the ALP couplings.
In addition to the contributions to the RG equations of the dimension-six SMEFT operators, the presence of the ALP also introduces modifications to the running of dimension-four couplings. Most of these modifications were discussed in [43]. However, additional contributions arise from the dimension-six operators included in (4), which we present here for the first time. They are where again Λ = 4πf .

Solving the ALP -SMEFT RG equations
The RG equations for the ALP -SMEFT Lagrangian can generically be written as 4 dC (4) where t ≡ ln µ, and we have collected the Wilson coefficients into vectors C (D) , with the superscript denoting the dimension of the associated operator, and the indices a, b (for D = 4), α, β (for D = 5), and i, j (for D = 6) labeling the corresponding coefficients. γ (D) are the corresponding anomalousdimension matrices, and γ (5,5) is the tensor accounting for the ALP source terms. The RG equations above form a system of coupled differential equations, for which an analytical solution is not known. One would thus be forced to solve this system numerically for a given set of initial conditions. It helps, however, to note that the RG equations above contain terms that are beyond dimension-six order in the EFT expansion. Indeed, using that C (4) (t) = C SM (t) + O(Λ −2 ), with C SM being the SM couplings, we can rewrite the RG equations for the Wilson coefficients of the higher-dimensional operators in the form in which all power-suppressed terms have been expanded out consistently. This system of equations admits an analytic solution for C (5,6) (t) in terms of C SM (t). 5 Explicitly, we find where t 0 and t f denote the logarithms of the initial and final energy scales, respectively, and the evolution tensors are defined as with T indicating that the exponential is t-ordered. Obtaining the evolution tensors is computationally expensive; however, they do not depend on the initial conditions set on the higher-dimensional operators and only need to be determined once (for a given set of SM input parameters). Once the evolution tensors are known, the computation of the running of the Wilson coefficients for arbitrary initial conditions is very fast. The evolution matrix U (6) for the SMEFT has already been determined and is part of the computer tool DsixTools 2.0 [72,73]. We have used a customized version of DsixTools, in which we have incorporated the ALP contributions to the RG equations, to compute the evolution tensor U (5,5) in the flavor-universal scenario. We provide this tensor for the SMEFT Wilson coefficients relevant for our fits in a Mathematica notebook as ancillary material. A new version of DsixTools featuring generic ALP contributions will be presented in a forthcoming paper.

Global analysis of ALP couplings
In the following, we utilize the ALP -SMEFT interference to derive (almost) mass-independent bounds on the ALP couplings to SM particles. Contrary to the direct bounds derived elsewhere, the results we obtain are model-independent in the framework we consider, consisting of an ALP added to the SM without additional sources of new physics. One of the advantages of this approach is that it lets us efficiently reuse results from existing SMEFT analyses.

Experimental inputs and SMEFT predictions
Limits on the SMEFT Wilson coefficients have been derived from a multitude of observables, including low-energy, Higgs and top data. For our global analysis, we utilize existing parametrizations of these observables in terms of the SMEFT Wilson coefficients and recast them in terms of ALP Wilson coefficients at the high-energy scale Λ. All of the SMEFT predictions used here are truncated at linear order in the dimension-six Wilson coefficients and employ {G F , α, M Z } as the electroweak input parameters.
Predictions for the Higgs sector are taken from [56] and references therein, while SMEFT predictions for the top sector are taken from fitmaker [55] and references therein. The experimental observables used for Higgs and top data are summarized in Tables 1-3 in Appendix A. The assumption of flavor universality in some SMEFT parametrizations of these observables, such as the ones from the experimental analysis in [74], causes certain complications, as ALP contributions are generally flavor-dependent. To overcome this issue, we will assume that the effects from operators involving quark couplings to gluons or the Higgs boson are dominated by those involving third-generation quarks, i.e.
where the notation [C x ] ij is used to denote the flavor indices i and j of the Wilson coefficient C x . The constraints on the remaining operators with flavor indices are typically experimentally dominated by first-and second-generation couplings. We therefore take a conservative approach and assume pure second-generation contributions for these operators in the considered Higgs and top observables, e.g. C (3) Hl ] 22 . We perform a χ 2 fit for the experimental data ⃗ d with the theory predictions ⃗ p(C i ) and covariance matrix V , using the definition For low-energy observables including electroweak precision data, neutrino scattering, atomic parity violation and quark pair-production at LEP2, we directly use the χ 2 function provided in [44,75], which includes the full flavor structure of the SMEFT Wilson coefficients. 6 The SMEFT predictions are translated to ALP predictions using the solution to the RG equations in (12) for C (6) i (t 0 ) = 0, i.e. neglecting possible matching contributions at the high scale. The value of the high-scale is set to µ 0 = Λ = 4πf , and the low-energy scale is identified with the relevant experimental scale for each of the observables, µ f = µ exp , with e.g. µ exp = m h for gluon-fusion Higgs production observables and µ exp = m h + 2m t for tth production. When considering ALP masses above the Z mass, as we will do in Section 3.3, we stop the ALP-induced running at m a and use the pure SMEFT running below this scale. ALP contributions to the RG evolution of dimension-four parameters such as the top-quark Yukawa coupling y t in (9) are found to be numerically irrelevant for the present analysis and have therefore been neglected.
The total χ 2 function is obtained by adding the individual contributions from low-energy observables, Higgs observables, and top data. The χ 2 functions for all these data sets, both written in terms of SMEFT Wilson coefficients and ALP parameters, are provided in the ancillary material.

Fit results
We present the limits on the ALP couplings at 95% CL and 99% CL in Figure 1. Assuming f = 1 TeV, we obtain O(1) bounds for C GG , C W W , C BB and C u , while C d and C e are much less

constrained, with limits of O(50).
Our choice of the scale f is motivated by the hope that new physics beyond the SM exists at a scale Λ = 4πf ≈ 10 TeV, as motivated by the hierarchy problem in light of current LHC results. Comparing the bounds from fitting one parameter at a time to those from a global analysis, we find that the limits on C W W , C BB , C d , and C e are minimally affected by the global analysis. Importantly, the weak constraints on C d and C e do not invalidate the limits imposed on the other Wilson coefficients in the global analysis. The global bounds on C u and C GG are weakened by 18% and 25% with respect to their one-parameter counterparts. For C GG , we notice that the corresponding limit favors a non-zero value at 95% CL, but it is compatible with zero at 99% CL. This discrepancy is caused by a minor experimental anomaly in three highly correlated bins of the CMS simplified template cross section analysis in the h → ZZ channel [76], which favors non-zero values of C uH or C uG in the SMEFT and consequently shifts C GG away from zero.
To investigate the correlations among the ALP couplings, we present two-dimensional fits in Figure 2. In each panel, we show the 95% CL bounds on two Wilson coefficients while setting the remaining coefficients to zero. As expected from the mild change of the global limits in Figure 1 with respect to the one-parameter fits, we find only weak correlations between most parameters. This is not surprising, since only few source terms of the SMEFT Wilson coefficients relevant for our analysis contain products of different ALP Wilson coefficients. Indeed, for the SMEFT Wilson coefficients appearing in our analysis, only C HW B as well as the dipole operators 33 contain the products C W W C BB , and C u C GG , C u C W W , C u C BB , respectively, in the leading-logarithmic (LL) approximation. The most interesting correlation patterns are observed for the combinations C GG -C u and C W W -C BB . For C GG -C u , a free-floating C GG allows C u to extend into a wider parameter region. The sign of the product of C GG and C u is however constrained to be negative at 95% CL when using the full data set. For C W W and C BB we find a slight preference for the product of the two coefficients to be positive.
It is interesting to study which of the considered data sets is the driving factor in constraining the individual Wilson coefficients. Individual bounds from low-energy, Higgs and top data sets are shown in Figure 2. Intriguingly, low-energy data dominate the constraints on C W W , C BB , C e , and even C u , which to first approximation describes the ALP-top coupling. As we discuss in the next subsection, the root of the strong constraint on C u from low-energy data is that it mixes under RG evolution with C HD , which is strongly constrained by the measurement of the W -boson mass. Only the bound on the ALP-gluon coupling C GG receives important contributions from Higgs and top data. There is an interesting interplay between the bounds from different experiments for C W W -C BB and C GG -C u . For C W W -C BB , Higgs data allow for a relatively wide parameter range as long as their product is positive. This degeneracy for C W W and C BB is broken by low-energy data.
For the pair C GG -C u , top data slightly favor a positive product of C GG and C u , while Higgs data (as well as the combination of Higgs and top data) favor a negative product.

Leading-log approximation vs resummation
It is interesting to investigate how the bounds obtained from an exact solution in (12), in which the leading-logarithmic corrections are resummed to all orders, differ from those derived using the LL approximation truncated at one-loop order, which leads to the simple formula For the low-energy observables at the Z pole, which dominate the fit for all coefficients except C GG , a full list of LL parametrizations can be found in Appendix B. We show in Figure 3 the oneparameter limits on the ALP couplings obtained using the one-loop truncated LL approximation alongside with the constraints obtained from a full (resummed) evaluation of the scale evolution. For C W W , C BB , C d , and C e , the LL approximation captures the dominant effects and the limits only change marginally with respect to the resummed evolution. However, for C GG and C u the LL approximation lacks important effects. To investigate this further, we display the limits on C u and C GG from different experimental sources in the two rightmost panels of Figure 3. We observe that the limits on C u from the top and Higgs sectors remain largely unchanged when using the LL approximation for the running. The main discrepancy arises from low-energy data, where the resummed running imposes tighter constraints than in the LL approximation. This discrepancy primarily originates from the RG evolution of the Wilson coefficient C HD . While the ALP contribution from C u to the Wilson coefficient C HD vanishes at LL order, it is simple to see that this is not the case for the resummed result. The full RG equation for C HD (neglecting for simplicity contributions proportional to α i ̸ = α s and all Yukawa couplings except for y t ) reads [67,68] where the ellipses refer to pure SMEFT contributions. Neglecting the running of the SM parameters results in the lowest-logarithmic approximation which is a two-loop effect. Thus, bounds on C HD only play a role in restricting the values of C u beyond the strict LL approximation. Even though this is a subleading effect, strong bounds on C HD from low-energy observables render it phenomenologically relevant for constraining C u . For C GG , shown in the right-most panel of Figure 3, the LL approximation is able to reproduce the constraints originating from the top sector quite well. Limits from low-energy data vanish completely in the LL approximation, but this data set only plays a marginal role in constraining C GG and hence does not influence the combined bounds. The dominant contribution to the combined constraints comes from Higgs data in this case. However, the LL approximation misses the most significant contributions from C HG and [C uG ] 33 for this data set, which are tightly constrained through gluon-fusion Higgs production and are only sourced by C GG beyond LL order. At lowestlogarithmic order, and taking into account the running of C u [64], the solutions for C HG , [C uG ] 33 in terms of C 2 GG are given by Both of the above estimates agree very well with the resumed results from the evolution tensor U (5,5) in (13).
Since most limits on the ALP couplings are well approximated at LL accuracy, we can deduce that they scale with the new-physics scale as |C i | ∼ Λ ln −1/2 (Λ/m Z ). As expected from the discussion above, the two exceptions to this scaling are C u , for which we find |C u | ∼ Λ ln −1 (Λ/m Z ), and C GG ∼ Λ ln n (Λ/m Z ) with n = −3/2 or n = −1, depending on whether [C uG ] 33 or C HG dominates.

Comparison with bounds from direct searches
In this subsection, we compare our indirect limits to direct limits on ALP couplings obtained in the literature. Since direct bounds are typically stronger for light ALPs, we focus on O(GeV) ALP masses, where we expect our indirect, (mostly) mass-independent limits to be more competitive. Direct limits in this regime are dominated by collider [31] and flavor [38] experiments. 7 We compare the direct and indirect bounds on the ALP Wilson coefficients in Figure 4, where the indirect bounds obtained from our global analysis are shown in red. The most relevant direct constraints on each of the individual ALP couplings are the following: • C GG : Flavor data impose the strongest direct constraints across the majority of the depicted ALP mass range. However, around m a ∼ 100 GeV, the inclusion of LHC multijet constraints becomes crucial [77]. Figure 4: Indirect 95% CL limits from ALP -SMEFT interference (red) compared to direct bounds from flavor, beam dump, and collider experiments, as well as supernova data. All direct bounds assume the remaining ALP Wilson coefficients to be zero. Direct bounds shaded in light gray are subject to additional model assumptions, see text for details. Note that we do not show the lower bound on C GG , excluding a non-zero value at 95% CL, as this bound disappears at 99% CL.
• C W W and C BB : In addition to flavor constraints, non-resonant ALP contributions to vector-boson scattering yield mass-independent direct bounds on ALPs with masses up to ∼ 100 GeV [78]. For non-resonant gluon-fusion ALP production followed by its decay to gauge bosons [79][80][81], constraints are placed on the products of the ALP-gluon coupling with the ALP-photon, C γγ = s 2 W C W W + c 2 W C BB , or the ALP-Z, C ZZ = c 2 W C W W + s 2 W C BB , couplings. 8 Specifically, |C ZZ C GG |/f 2 < 4 · 10 −2 TeV −2 [81] and |C γγ C GG |/f 2 < 5 · 10 −3 TeV −2 [79]. As these bounds involve products of ALP parameters, they are not shown in our plots which depict the bounds for one coupling at a time. However, we point out that they provide the dominant bounds in the limit in which both Wilson coefficients are sizable. In the high ALP-mass region, we present bounds derived from collider constraints on the ALPphoton coupling [32]. The corresponding limits are shown in light gray to highlight their model dependence. It is important to note that for heavy ALPs with m a > m Z , where additional decay channels like a → Zγ open up, a dominant decay to photons becomes highly unlikely [82]. In darker gray, we present the same limits assuming a branching ratio into photons of 10 −3 .
• C u and C d : Constraints on the ALP parameter space stem from flavor data that, for the case of C u , cover the full displayed ALP-mass range. In addition, C u gets further constrained by LHC tt searches for m a ≤ 100 GeV [83].
• C e : We present direct bounds from flavor and dark photon searches at BaBar [84], which cover a similar mass range compared to constraints from Υ decays. Furthermore, we include constraints from SN1987A supernova observations [85] and beam dump searches at SLAC [86], which are relevant in the m a < 1 GeV mass range. In the 11.5 − 50 GeV mass range, we consider LHC constraints on h → aµμ [87], assuming a 100% branching ratio of the ALP to muons. These constraints are shaded in a lighter gray to account for the possibility of decays to taus, which would weaken the bounds.
Overall, we find that the ALP -SMEFT interference can constraint previously untested regions of the ALP parameter space. Furthermore, most direct bounds have specific model assumptions, often requiring all remaining coefficients to be zero, that do not apply to our indirect bounds. The indirect limits presented here thus offer good complementary probes, even in cases where the direct limits would a priori seem more competitive.

Interpretation in terms of UV-complete models
We now reinterpret our bounds in terms of UV-complete ALP models. In particular, we describe the effect of including SMEFT contributions beyond the RG-induced effects considered before, which stem from the incorporation of (tree-level) threshold corrections. We focus in the two main (fundamental) axion UV completions: the Kim-Shifman-Vainshtein-Zakharov (KSVZ) [39,40] and the Dine-Fischler-Srednicki-Zhitnitsky (DFSZ) [41,42] models, which have been already proposed as ALP benchmark scenarios [88].

KSVZ model
The KSVZ model extends the SM with a pair of fermions, Q L and Q R , which transform non-trivially under SU (3) c and chirally under a global U (1) A symmetry, and a SM-singlet scalar S, which is only charged under the U (1) A symmetry and acquires a non-zero vacuum expectation value (vev), thus spontaneously breaking the global symmetry. The most general Lagrangian for this model reads where y Q , µ S , λ S , and λ SH are real parameters and L Qq is a possible portal coupling between Q and a SM fermion. As discussed in [89], this term is introduced to let the extra quarks decay, as otherwise the Lagrangian would be invariant under a vectorial U (1) Q symmetry that would make them stable. In the original KSVZ implementation, where the extra quarks transform under the SM as Q L,R ∼ (3, 1) 0 , no renormalizable term for L Qq is possible. However, there are multiple representation choices for which L Qq ̸ = 0. For concreteness, we consider here the case where Q L,R ∼ (3, 1) −1/3 and the global U (1) A charges are X S = 1, X Q L = 1, and X Q R = 0. In this case, we have Additionally, we consider possible soft U (1) A -breaking terms that will give mass to the pseudo-Goldstone ALP at the expense of spoiling the solution to the strong CP problem. For definiteness, we consider the term with κ being a real parameter controlling the size of the breaking. In the U (1) A -broken phase, it is convenient to write the scalar singlet as with f denoting the vev of S, and ρ and a corresponding to the radial and pseudo-Goldstone components of the field, respectively. After performing the fermion shift which removes the ALP from the Yukawa interaction, the UV Lagrangian in the U (1) A -broken phase reads with M Q = y Q f / √ 2, M 2 ρ = λ S f 2 and m 2 a = 2κ 2 . Note that, after U (1) A -symmetry breaking, the Higgs doublet mass gets a correction of order λ SH f . Thus, if the hierarchy between the electroweak scale and f is large, one would expect λ SH ∼ v 2 /f 2 to avoid a fine-tuned cancellation of parameters.
We consider the scenario in which M ρ , M Q ∼ f are heavy and integrate out the corresponding particles. The resulting EFT Lagrangian at tree-level order and up to dimension-six interactions reads 9 where the second and third terms are removed by an appropriate redefinition of the Higgs potential parameters. Namely, Furthermore, we see that the explicit U (1) A -breaking term not only gives mass to the ALP, but also introduces other shift-symmetry breaking interactions, m 2 a a 4 and m 2 a a 2 (H † H). 10 Finally, Q i are dimension-six SMEFT operators whose definition can be found in [63].
We now turn to analyzing constraints on the KSVZ model from Higgs, top and low-energy data. As discussed above, the KSVZ model features a fermiophobic axion (at tree-level order), with different Q charges under the SM gauge group yielding different values C GG , C W W and C BB . A common feature of all KSVZ models is the presence of a non-zero Q H□ , which is generated via the tree-level exchange of the associated scalar radial excitation. Profiling over the remaining parameters, we obtain the constraint λ 2 S f /λ SH > 2.8 TeV from the limits on the Q H□ coefficient. On the other hand, the limits on the bosonic ALP couplings when profiling over the remaining parameters (including the Wilson coefficient of Q H□ ) do not change by more than 10% with respect to the one-parameter limits presented in Figure 1. Therefore, we refer to this plot for limits on KSVZ models with generic Q charges.
For KSVZ models with additional portal couplings for the heavy vector-like quarks, such as the one presented in (27), we can additionally set constraints on the coupling strength of the portal coupling y q . Assuming for simplicity that y q is flavor universal, we find the limit |y q /M Q | < 0.1 TeV −1 , which is dominated by the strong constraints on C

DFSZ model
The DFSZ models consists of a two-Higgs-doublet, H 1,2 , plus SM-singlet, S, scalar extension of the SM. The Lagrangian of the model is chosen such that it preserves, at the classical level, a global U (1) A symmetry and reads where ⊃ denotes that we omitted the SM-like terms in the Lagrangian. All scalar potential parameters in the Lagrangian above are real, including λ SH 12 which can be made real by appropriate global phase redefinitions of the fields. In the charged-lepton Yukawa, i = 1, 2 corresponds to two different versions of the model, which we denote as DFSZ I and II, respectively. The last term also admits a different choice, with S rather than S 2 , corresponding to a different U (1) A charge implementation. Different choices for this term have mild effects in the ensuing discussion and we focus on this variant of the model for definiteness. As we did for the KSVZ model, we admit the possibility of an explicit U (1) A -breaking term, which we choose to be the same as before with κ being a real parameter controlling the size of the breaking. As before, the scalar potential parameters are chosen such that S acquires a vev that breaks the global U (1) A symmetry. Once more, we parameterize the SM-singlet as The two-Higgs-doublet spectrum can be rather different depending on the values of m 1 , m 2 and λ SH i . Here, we assume that these parameters are such that we are in a decoupling regime where a full doublet and ρ are much heavier than the ALP and the SM particles. 11 The heavy doublet, Φ, and SM Higgs, H, are linear combinations of H 1 and H 2 . Namely, where m 2 ii = m 2 i + λ SH i f 2 /2 (i = 1, 2), m 2 12 = λ SH 12 f 2 /2, and R(α) is a 2 × 2 rotation matrix. Once this rotation is taken into account, the SM Yukawas are related to the mixing angle and the original Yukawas as with c α ≡ cos α, and s α ≡ sin α. Perturbativity constraints on the UV Yukawa couplings restrict the possible values of α. Using the perturbative unitarity bound |Γ 33 u | ≲ 3 [93], we get the constraints |c α | ≳ y t /3 and |s α | ≳ y b /3, independently of the DFSZ type.
The pseudo-Goldstone, a, can be moved away from the scalar potential, yielding a coupling structure like the one in (2), by means of the following shifts of the scalar and fermion fields with X H 1 = 2s 2 α and X H 2 = −2c 2 α . 12 After performing these shifts and the Higgs rotation in (32), we obtain the following UV Lagrangian in the U (1) A -broken phase where M 2 ρ = λ S f 2 , m 2 a = 2κ 2 , t α ≡ tan α and η α = −t α (t −1 α ) for DFSZ I (II), and we omitted the Lagrangian terms involving heavy fields (either Φ or ρ) that do not contribute to the tree-level matching at dimension six. The dimension-five ALP couplings depend on the model variant and are given by Finally, we have defined the following couplings for simplicity We integrate out the ρ and Φ fields at tree-level with the help of Matchete [90]. The resulting EFT Lagrangian at dimension-six reads This choice of H1,2 shifts is uniquely determined by the requirement of having no contributions to the ∂ µ a (H † i ← → DµH) operator, which would introduce a mixing between the ALP and the Z would-be Goldstone boson after electroweak symmetry breaking. Figure 5: Two-dimensional limits on α vs f in the DFSZ model for two benchmark scenarios: S1 (benchmark) and S2 (profiled), see text for details. The dark and light gray bands corresponds to regions where pertubative unitarity is violated (namely, Γ 33 u ≳ 3) and where Γ 33 u ≳ 1, respectively.
where C ψH and C H are given in terms of the original scalar-potential parameters by Analogously to what we did in the KSVZ case, we have redefined the SM Higgs potential parameters, µ and λ, to account for the matching corrections. If the hierarchy between the electroweak scale and f is large, one would again expect λ SH ∼ v 2 /f 2 to avoid a fine-tuned cancellation of scalar-potential parameters.
Translating the above notation to the one used in our global analysis (cf. (4)), we find that only ALP couplings to fermions are non-zero in both DFSZ models: C u = −C e = −2s 2 α and C d = −2c 2 α for DFSZ I, and C u = −2s 2 α and C d = C e = −2c 2 α for DSFZ II. We present the bounds on the mixing angle α and the ALP decay constant f in Figure 5. We consider two scenarios: S1) where the coefficients of the SMEFT operators Q H , Q H□ , and Q ψH are assumed to be suppressed, as would be expected if the scalar-potential parameters are small; and S2) where we profile over the coefficients of these operators within the range |C ψH |, |C H□ |, |C H | < 1. Since the ALP couplings to fermions in the DFSZ I and DFSZ II models differ only by their couplings to leptons, which are weakly constrained, we find that the limits on both models are (almost) identical. As shown in Figure 5, the overall limit on f is dominated by the matching correction from the four-quark operators Q (8) qu and in particular Q (1) qu , which run into the coefficients of the SMEFT operators Q HD , Q (1) Hq and Q (3) Hq that are tightly constrained at the electroweak scale. The obtained limits are found to be weak, except in the region where the UV Yukawa couplings are larger than one and dominates close to the non-perturbative Yukawa region. Limits from the ALP coupling to up-type quarks are suppressed by s 2 α and thus only play a subdominant role. When profiling over the other matching corrections within the range |C ψH |, |C H□ |, |C H | < 1, we observe that the limits on both DFSZ models are slightly diminished with respect to S1, especially for intermediate values of |α|.

Conclusions
While the SMEFT is normally used to describe the possible effects of yet undiscovered heavy particles, we have shown in this paper that SMEFT analyses can also be reinterpreted to infer indirect information on light new physics. Exploiting the non-trivial RG flow of the ALP couplings into the SMEFT Wilson coefficients [43], we have used existing SMEFT studies including lowenergy, Higgs and top data to constrain these couplings. Contrary to other ALP constraints, the ones presented here posses the unique feature of being largely independent of particular assumptions on the ALP mass, lifetime or branching rations.
Furthermore, we have obtained for the first time, a semi-analytic solution to the ALP RG equations at dimension six under the assumption of flavor-universal ALP interactions. This solution, given in the form of RG evolution tensors, can readily be used to derive ALP and SMEFT couplings at an arbitrary scale (provided there are no mass thresholds). Its generalization to generic ALP interactions and the incorporation of threshold corrections into a modified version of DsixTools will be presented in a forthcoming publication, thus paving the way for automated ALP analyses. Even with the present assumptions, the ALP-to-SMEFT evolution tensor provided in the ancillary files can readily be used in most SMEFT analyses, such as the one presented in this paper.
The resulting bounds on the bosonic ALP interactions C GG , C W W , C BB and the ALP coupling to up-type quarks C u are found to be of O(1) for f = 1 TeV. The couplings to down-type quarks and leptons, C d and C e , remain weakly constrained, with limits of O(50), as expected given the additional Yukawa suppression present in these couplings. The bounds on C W W , C BB , C u , C d , and C e primarily arise from low-energy precision observables, such as measurements at the Z pole. On the contrary, the limit on C GG is mainly driven by Higgs and top physics. In our global analysis, we found weak correlations between Wilson coefficients. However, there is an interesting interplay between the limits on C u and C GG , where a negative product of the two Wilson coefficients is favored, and on C W W and C BB , where any large product of the two coefficients is ruled out. We have also found that the LL approximation captures most of the phenomenologically-relevant effects for all ALP couplings except C u and C GG , which generate non-trivial contributions to strongly constrained SMEFT directions at higher order in the RG resummation.
Comparing with direct ALP searches, our limits constrain large regions of previously uncovered areas of the parameter space for ALP masses above 10 GeV. Even for lower masses, the obtained bounds can partly compete with existing ones and have the major advantage of being independent of specific assumptions on the ALP properties, thus offering a complementary probe in regions where direct bounds would a priori dominate. It would be interesting to investigate how the combination of direct and indirect bounds further narrows the ALP parameter space in a global analysis. We leave the comprehensive study of both direct and indirect constraints in a global analysis for future work.
We have also reinterpreted our results in the context of two benchmark UV completions based on the KSVZ and DSFZ models. The KSVZ model features no tree-level couplings to fermions and we have found that tree-level threshold corrections, arising from integrating out additional heavy particles present in the model, do not significantly affect the limits on C GG , C W W and C BB obtained from the ALP -SMEFT analysis. On the contrary, the DFSZ models we studied yield only fermionic ALP couplings. In this case, we found that tree-level threshold corrections can be more important and dominate the model constraints. However, the model remains weakly constrained except in regions where the UV Yukawas are large. A more dedicated study including different assumptions on the mass spectrum, additional observables or incorporating one-loop threshold corrections remains to be explored in future studies.
Hq ij = − Hq ij = − Hl ij = − where again Λ = 4πf , and the SMEFT Wilson coefficients normalized as C i /Λ 2 . Parametrizing the Z and W coupling modifications as in [44], we obtain the following relations: with the additional relations δg Zν L = δg Ze L + δg W l L and δg W q L ≈ δg Zu L − δg Zd L (for V CKM ≈ 1). From these expressions it becomes clear why C u remains unconstrained by Z-pole observables when working at LL accuracy, as this parameter only enters in Ztt coupling modifications.
[102] CMS collaboration, Combined Higgs boson production and decay measurements with up to 137 fb-1 of proton-proton collision data at sqrts = 13 TeV, CMS-PAS-HIG-19-005 (2020) . [135] ATLAS collaboration, Inclusive and differential measurement of the charge asymmetry in tt events at 13 TeV with the ATLAS detector, . [142] CMS collaboration, Measurement of differential tt production cross sections in the full kinematic range using lepton+jets events from pp collisions at √ s = 13 TeV, .