Phenomenology of axion-like particles with universal fermion couplings – revisited

Axion-like particles (ALPs) emerge in many extensions of the Standard Model as pseudo-Goldstone bosons of a spontaneously broken global symmetry. Understanding their phenomenology in high-energy collisions is crucial for optimizing experimental searches and understanding the exploration potential of future experiments. In this paper, we revise the phenomenology of ALPs with universal couplings to fermions. In particular, we analyze the hierarchy and uncertainty of the various ALP production channels depending on the proton collision energy and the placement of the experiment, and provide improved calculations of the hadronic decay modes.


I. INTRODUCTION
Axion-like particles (ALPs) a are pseudoscalar particles that arise in theories with spontaneously broken global chiral symmetries, generalizing the idea of the QCD axion -a hypothetical light particle capable of solving the strong CP problem [1][2][3].While the QCD axion obtains its mass directly from its coupling to gluons, a generic ALP may interact with various particles and have an arbitrary mass [4].The lowest-order gauge-invariant Lagrangian of the ALP interaction takes the form [5][6][7] where G, W, B are field strengths of the Standard Model SU c (3), SU L (2), and U Y (1) gauge groups, α s , α W , α B = g 2 /4π are corresponding running couplings, Gµν = 1 2 ϵ µναβ G αβ is the dual strength and Ψ F denote the SM fermion multiplets.Furthermore, F is a dimensional scale, C G,W,B are dimensionless parameters, and C F are hermitian matrices characterizing the structure of the ALP couplings to fermions.For light ALPs with m a ≃ 1 GeV, past experiments have excluded combinations F/C i ≲ 1 TeV [8,9].Therefore, ALPs belong to the class of the so-called Feebly-Interacting Particles, or FIPs.Many studies have explored the phenomenology of the individual terms in the Lagrangian above, in particular couplings to gluons [10][11][12], photons [13][14][15], electroweak gauge bosons [16][17][18], fermions [19][20][21], the effect of flavour-violation [22,23] and renormalisation group evolution [24][25][26][27], as well as the interplay between different couplings [28][29][30][31].One case of particular interest is ALPs that interact dominantly with fermions with universal and flavourdiagonal couplings, where ℓ are leptons and q are quarks [32,33].The case C ℓ = C q has been identified by the Physics Beyond Colliders (PBC) initiative [8] as one of the benchmark models (called BC10) to demonstrate the FIP exploration potential of future experiments.Due to a lack of in-depth theoretical studies of this model, the description of the phenomenology of such ALPs proposed in Ref. [8] suffers from two issues.First, following Ref.[19], the contribution of the hadronic decays to the total ALP decay width is assumed to be negligible.While this is a reasonable assumption for light ALPs with m a ≪ 1 GeV, where the only relevant decay mode is into three pions [25], the hadronic decays may actually dominate the total width for heavier ALPs, as indicated by calculations of the ALP decay width into quark and gluon pairs [25,27].Another problem is that the production of such ALPs at beam dump experiments has been approximated by considering decays B → X s + a, where X s = K, K * only [20,27].Nevertheless, there may be a sizable contribution from the B decays into other resonances, X s = K 1 , K * 2 , K * 0 .For the case of light Higgs-like scalars [34], which couple to the b → s operator in a way similar to ALPs, these decays have been shown to correspond to 1/3 of all possible decays.The same effect can be expected to be relevant for the case of ALPs with interactions as in eqs.(1) and (3).Second, additional production processes arise due to the mixing of the ALPs with light pseudoscalar mesons m 0 = π 0 /η/η ′ , similar to the ALPs coupled to gluons.
In this paper, we address these issues by revising the phenomenology of the model given in eq. ( 3).Our goal is to provide a detailed and comprehensive description of the PBC BC10, which the community may easily implement to consistently interpret ALP constraints from existing experiments and study the projected sensitivities of proposed searches.The results summarized in this paper are also accessible in a Mathematica notebook supplementing the paper (see Appendix B). 1 We implement the revised phenomenology in SensCalc [35] -a public code that uses the semi-analytic approach to calculate the number of events with decays of FIPs and sensitivities of proposed experiments.
The remainder of this work is organized as follows.In Sec.II, we discuss the details of the ALP model that we consider, and in particular, the choice of the scale at which the underlying Lagrangian is defined.In Sec.III, we discuss various contributions to the ALP production flux depending on the proton collision energy for the given experiment and its geometric placement.In Sec.IV, we study the decay palette of the ALP.We conclude in Sec.V.

II. MODEL DETAILS
The phenomenology of GeV-scale ALPs depends on the scale Λ at which the interactions in eqs.(1), (2) are defined.This is because of a non-trivial renormalization group (RG) flow, which generates additional effective couplings absent in eqs.(1), (2) and modifies the initial couplings C i if the energy scale µ of the process under consideration differs from Λ.
Typically, Λ and F are closely related, F ≃ Λ [25].The values of the Wilson coefficients C i (Λ) at that scale, however, are in general unknown, so the quantities f i ≡ F/C i (Λ), controlling the interaction strength of ALPs with the corresponding SM fields, can be treated as independent parameters.Here and below, we will introduce the coefficients c i (µ) ≡ C i (µ)/C i (Λ), such that by construction c i (Λ) ≡ 1, and perform our studies in terms of the effective suppression scales f i , assuming that f i > 0 for definiteness.The assumption of universal couplings to all the fermions then corresponds to a single scale f f = f .With this assumption and setting Λ = 1 TeV, the model (2) exactly coincides with the BC10 model from [8], such that our results may be directly used by 1 Available on https://github.com/maksymovchynnikov/ALPsphenomenologythe experimental community: For a specific ultraviolet model that predicts the Wilson coefficients C f (Λ), our results can be directly reinterpreted in terms of the fundamental scale F .The RG evolution has been thoroughly studied for general models of ALPs (see Refs. [25][26][27] and references therein).The evolution is usually split onto the flow from Λ down to the electroweak scale µ w ≡ m t and from µ w down to the scale of the process with ALPs, which is of the order of the ALP mass.For the processes with hadrons, there is one additional scale ∼ 4πΛ QCD where the perturbative QCD should be matched with Chiral Perturbation Theory (ChPT).
The couplings C G , C W , C B entering the Lagrangian in eq. ( 1) are invariant under the RG evolution at least up to the second order in the loop expansion.Therefore, they will only be generated from the initial Lagrangian in eq. ( 3) through threshold corrections when integrating out heavy quarks.This is not the case for the fermion couplings from eq. ( 3), which may also evolve due to electroweak and strong interaction loops.To study the dynamics of these couplings, we solve the RG equations from Ref. [25].As for the couplings to heavy quarks q = c, b, t, we define that they do not run below the scale µ = m q as the corresponding degrees of freedom are integrated out; in particular, c t (µ < µ w ) = c t (µ w ).In Fig. 1 (left panel), we show the value of the running couplings c i (µ) from eq. ( 3) to various SM fermions at the scale µ = 2 GeV as a function of Λ.
The RG dynamics of the lepton couplings c ℓ is driven by electroweak interactions.The ratio (c ℓ (Λ) − c ℓ (µ))/c ℓ (Λ) is ≲ 0.2 for Λ ≲ O(10 TeV).Since all the leptons have the same properties under the electroweak gauge group (such as the weak isospin and electric charge), the flow of c ℓ is lepton-flavor independent.This is not true for the quark couplings.Their evolution down to µ w is flavor universality violating because different quarks have different weak isospins and electric charges.Assuming Λ = 1 TeV, the relative difference between u and d, s couplings is (c d/s (µ) − c u (µ))/c d/s (µ) ≃ 10%.As we will see in Sec.III, this difference cannot be neglected when studying the interactions of the ALPs with neutral pions (see Sec. III).
The coupling universality between quarks and leptons gets broken even if assuming Λ = µ w .This is because below µ w the flow of the c q s is determined by loops involving strong interactions, resulting in the corrections of the order of 10 −2 c q , whereas c ℓ s evolve only due to the EM corrections, which are of the order of 10 −5 c ℓ [26].
The ALP-gluon interaction is another type of interaction important for ALP production and decay.At the leading order in α s , the matrix element of the type GG → a (the gluon fusion process) or a → GG (the ALP decay into a pair of gluons) is generated by the following matrix element: where G a µν (p) is the linear part of the gluon strength tensor with the replacements ∂ µ → ip µ and G a µ → ϵ a µ (p), with ϵ a µ being the polarization vector.The effective cou- with and τ f = 4m2 f /m 2 a .In the regime τ ≫ 1, the function B 1 behaves as B 1 (τ ≫ 1) ≈ −(3τ ) −1 .In the opposite regime, B 1 (τ ≪ 1) ≈ 1.Therefore, the coupling to gluons is mainly generated by the quarks lighter than the ALP.
Additionally, the RG flow (via loops involving top quarks and charged weak gauge bosons) generates the flavor-changing neutral current (FCNC) coupling D i D j a, where D = d, s, b are down quarks and i ̸ = j [26]: where ) and c ij is the model-dependent coupling: with The term I t (Λ, µ w ) represents the contribution of the RG flow from the scale Λ down to µ w = m t and vanishes if Λ = µ w .
Neglecting the mass of the lighter quark in eq. ( 7), the FCNC Lagrangian may be rewritten as [36] where , where m qj is the mass of the heavier quark among the pair q i q j .The FCNC couplings for the transitions s → d and b → s have a huge impact on the ALP phenomenology as they may dominate the production of the ALPs depending on the amounts of kaons and B mesons produced in the given experiment.The value of this coupling is very sensitive to Λ, growing by two orders of magnitude if increasing the scale from Λ = µ w to Λ ∼ 10 TeV (see Fig. 1, right panel). 2  Taking this into account, we consider two representative choices of Λ: the one with Λ ≡ µ w , and another one with Λ ≡ 1 TeV, which is the reference scale used for the PBC BC10 benchmark [8].

III. ALP PRODUCTION
The ALPs in eq. ( 3) may be produced by decays of kaons and B mesons, via the mixing with light pseudoscalar mesons m 0 , or by deep inelastic scatterings (DIS).
We describe these production channels in detail in the section below.The amounts of the produced mesons and the DIS cross-section are collision energy dependent.Therefore, to make an experiment-independent   9) assuming the model (3) with f = 1 PeV and two scales at which it is defined: Λ = mt, and Λ = 1 TeV.
comparison, we will consider three collision energies √ s ≈ 16 GeV, 28 GeV, and 13 TeV (corresponding to the collision energies at DUNE, the SPS beam and the LHC).We took the meson production fractions from the SensCalc repository [35].

A. Decays of B, K mesons
We will consider the interactions abs, abd, and asd, for which the quark running in the loop is the top quark; the other interactions are heavily suppressed by the Yukawas of lighter quarks and/or CKM elements and are irrelevant.The corresponding decay processes are B → a+X s , B → a + π, and K → a + π.As we see from eq. ( 8), the values of the couplings describing these transitions differ only by the CKM products V * ti V tj .This product is the largest for the b → s transition; however, the other processes are also important.Namely, the abd coupling is suppressed by |V t→d /V t→s | ≈ 0.2; however, the process B → a + π is the only possible above the threshold B → K + a.The relative suppression of the sd coupling is even larger, |V t→d | ≈ 8 • 10 −3 .However, depending on the experiment, the number of kaons may be much larger than that of B mesons, which may compensate for this suppression.The values of the corresponding couplings for the two different scales Λ = 1 TeV or Λ = µ w are given in Table I.In particular, the value C bs (1 TeV) matches with the one used for the BC10 model [8,36].
Having the operator of the FCNC interaction (9), one may calculate the matrix elements of the processes B/K → a + X, where X is a hadronic state containing an s quark or a d quark.They have the form where the parity-even and parity-odd transition matrix elements are Because of the parity conservation in QCD, if m, m ′ have the same parity, only M (S) m→m ′ contributes, while for m ′ having a different parity than m only the The matrix elements (11) match with the matrix elements M XX ′ from eq. (B.7) from [34], used to compute the production of the Higgs-like scalars, which is caused Br(B->a+X) where X is a hadron that contains an s or d quark.By the K * 0 channel, we denote the final states K * 0 (700) and K * 0 (1430), by K1 -K1(1270), K1(1400).The dashed black lines correspond to the probability of the process B → K/K * (892) + a considered previously in the literature.
by the similarity of the FCNC operator for ALPs and Higgs-like scalars [37][38][39][40].Therefore, instead of computing the branching ratios using eq.( 10) one may use the results of Ref. [34] after the rescaling of the branching ratio with the proper coupling.
Ref. [34] used the matrix elements computed using light-cone QCD sum rules and considered the mesons 2 , and X d = π.The branching ratios of various decays are shown in Fig. 2. Compared to the literature where only the decays B → a + K/K * (892) have been considered [8,10,29], we find almost 4 times larger total production probability.In particular, in the domain of masses m a ≲ 1 GeV, the dominant production channel is into K 1 and K 0 * mesons.

Interaction
If the ALP is light (m a ≲ 2 GeV), the description of its hadronic interactions in terms of the qq and GG operators from eq. ( 1) becomes inadequate since the QCD enters the non-perturbative regime.Instead, light mesons and their interactions represent the strongly interacting sector.
We follow the existing studies [10,25,26] and obtain the Lagrangian of the ALP interactions with the pseudoscalar mesons P = π, η, K, . . .by using the matching of the operator (3) with the ChPT Lagrangian.
Details are provided in Appendix A; we summarize the main features below.In general, the interaction (3) leads to the kinetic mixings of the ALP with neutral pseudoscalar mesons m 0 .This contrasts with the case of the ALPs coupled to gluons, where the gluon operator also induces the mass mixings [10].We need to diagonalize the kinetic term to find the relevant interactions.The fields m 0 entering the Lagrangian are related to the mass eigenstates m 0 phys , a phys by Here, the second term is due to the kinetic mixing with the ALPs, and the third one appears from the mass mixing between the mesons emerging from the minimal ChPT breaking term.
In the limit m s ≫ m u,d , the mixing angles are is the isospin symmetry breaking parameter, f π ≈ 93 MeV is the pion decay constant, and F (m a ) is a phenomenological function ensuring the drop of the VMD contribution according to quark counting sum rules [10,29].Similar to [10,41], in eq. ( 13), we fix θ η = arcsin(−1/3), motivated by the fact that various phenomenological studies of the effective Lagrangian of the decays of light mesons consider this value.The expressions (13) are given in terms of the couplings c u , c d , c s instead of a single c q to account for the RG flow (remind Sec.II).
In the first order on the parameter f π /f ≪ 1 and far from resonance domains m a = m m 0 , the ALP and meson masses are left unchanged.Therefore, the only impact of the diagonalization (12) is the appearance of new interactions between the ALPs and mesons.
The ALP mixing with π 0 emerges either from the RG flow (the first term in eq. ( 13)), or via the mixing of unphysical π 0 with η and η ′ (the second term).
The behavior of the squared mixing angles as a function of the ALP mass for the two reference scales Λ = m t and Λ = 1 TeV is shown in Fig. 3.

Production
All relevant production processes with ALPs and light neutral mesons occur through their mixings (13).We assume that the ALP production cross-section due to the mixing is given by However, depending on the ALP mass, its kinematics would be different from the corresponding meson kine-matics.We follow the procedure described in [29] to account for this.
We should stress that this description does not account for the ALP mass dependence of the production crosssection.Clearly, the production of the ALPs heavier than m 0 at the unit mixing angle must be kinematically suppressed, which is not considered in eq. ( 14).This point is to be improved in future works.

C. Deep inelastic scattering process
Another important process of ALP production is deep inelastic scattering.At the parton level and at the leading order in α s , it is described by the fusion where for the second process, the matrix element is given by eq. ( 4).The parton model applicability breaks down if the characteristic scale of the process √ s qq = m a becomes comparable with Λ QCD .We "turn on" the process (15) at m a = 1.5 GeV.The DIS process is the only relevant production channel for heavy ALPs with m a > m B − m K , given that its kinematic threshold is extended until the center-of-mass energy at the experiment.
To estimate the cross-section of the process (15), we implement the fermionic and gluonic matrix elements in MadGraph5 [42] using FeynRules [43,44] and generate the leading-order processes of the gluon and quark fusion.For the parton distribution function, we use NNPDF 3.1 NNLO set, which is a common choice for FIP sensitivity studies [45].
We have found that the quark fusion is strongly suppressed compared to the gluon fusion.The reason for this is a large value c eff G ≈ c u + c d + c s ≈ 3 and the fact that the gluonic squared matrix element is proportional to m3 a rather than to m a m 2 q as in the case of the quark fusion, where m q ≪ m a .
The cross-section of the gluon fusion at various centerof-mass collision energies is shown in Fig. 4. It suffers from a large systematical uncertainty, which we illustrate by varying the renormalization and factorization scales µ r = µ f = m a by a factor of 2. A huge fraction of this uncertainty comes from the scaling σ gg→a ∝ α 2 s , which varies rapidly for small scales µ = O(1 GeV).Another important question relevant for the LHC is that the production of light ALPs sits at a very small energy fraction x = m 2 a /s pp , where the PDFs have a large uncertainty.Including the uncertainties from PDFs would further increase the overall error: by considering several other choices of PDFs, we have found that the crosssection differs by an additional O(50%).
Unlike the production via mixing and decays of B mesons, the DIS cross-section is practically independent of the RG flow.Indeed, the effective gluon coupling includes the summation over u, d, s quarks.For the ALP masses of interest, it is proportional to the combination (c u + c d + c s ) 2 , which for the wide range of the choices of the scale Λ changes insignificantly.

D. Discussion
To compare the contribution of the particular production channels to the ALP yield, we consider the production probabilities per proton collision: Here, χ X is the fraction of the produced X particles per proton collision (for B, we take both mesons and anti-mesons), and σ pp is the total proton collision crosssection.The mixing angles θ m 0 a are taken from eqs. (13), and for the branching ratio Br(B → X s a) we used Fig. 2 and Table I.We will not consider here the ALP production by decays of kaons, since the latter are long-lived, and the ALP flux may be heavily affected by the interaction of the kaons with the infrastructure surrounding the kaon production point.The most important factor is their absorption by the material: the absorption length of the kaons is typically smaller than their decay length 3.7E K /m K m.For the beam dump experiments with thick targets, kaons would already be heavily absorbed inside the target.For the LHC-based experiments, the effective kaon decay volume is limited by the detector, and only a small fraction of kaons would decay there due to their huge boosts (see, e.g., [46]). 4he total production probabilities for these energies are shown in Fig. 5. From the comparison, we see that the dominant production channels at the LHC are decays of B mesons (thanks to the large fraction of the produced b b pairs) and Drell-Yan process, independently of the model scale Λ chosen.This finding differs strongly from the case of ALPs coupled to gluons [12,25], for which the  16)) per proton collision at the LHC, SPS, and FNAL (DUNE), with the collision energies of √ s = 13 TeV, ≈ 28 GeV and ≈ 16 GeV correspondingly, for the model (3).The probabilities are evaluated for the value of the ALP coupling f = 1 PeV and assuming two scales Λ: Λ = mt (dashed lines) and Λ = 1 TeV (solid).The meaning of the vertical gray bands is the same as in Fig. 3. See text for details.
sensitivities to FIPs, therefore, requires a dedicated study (see, e.g., Ref. [49].main production is via the mixing with neutral mesons.The reason may be understood in the following way.All the production channels except for the flavor-violating meson decays receive similar contributions from the treelevel couplings to fermions and gluons c f , c G , while the rates of the FCNC decays are very different.Namely, the main contribution to the FCNC coupling is made via the ALP coupling to the top quark, which is absent at the tree level for ALPs coupled to gluons.Instead, it is generated by loops involving gluons [12].The ratio between bs couplings in cases of the fermionic and gluonic ALPs is of the order of f G /f (α s /π) 2 ln(Λ/m t ), which is ≫ 1 if assuming similar couplings f and f G ≡ F/C G .
The mixing with light mesons may become relevant for the experiments operating at lower energies.Depending on the scale Λ, it may dominate the total production at SPS at masses m a ≲ 2 GeV.It also dominates the production at FNAL, even at large Λ, given the small center-of-mass energy and the correspondingly tiny fraction of produced B mesons.
To conclude this discussion, we emphasize that the hierarchy of the production channels may change depending on the placement of the experiment with respect to the proton beam axis.Generically, the angular distributions of the light ALPs from B decays and those produced by the Drell-Yan process are broader than the distribution from the mixing with light mesons.Therefore, if these production channels provide similar overall amounts of ALPs, the mixing with the mesons would dominate for on-axis experiments with small angular coverage, while the other channels dominate for off-axis experiments.To demonstrate this point, in Fig. 6, we show the production probabilities for the ALPs flying within the polar range θ < 10 mrad and θ > 10 mrad at SPS, assuming Λ = 300 GeV.

A. Decays into leptons and photons
The matrix element of the decay of ALPs into a pair of photons is given by Depending on the ALP mass, the effective coupling where m ′ a ≃ 2 GeV is similar to the matching scale between the ALP's ChPT and QCD perturbative decays (see the next subsection).
The VMD contribution c VMD γ originates from the mixing of the vector mesons ρ, ω, ϕ with photons.The corresponding contributions are We have checked that this coupling approximately reproduces the widths of the anomalous decays m 0 → 2γ with m 0 = π 0 /η/η ′ in the symbolic limit The loop contribution is given by triangle diagrams with fermions f = l, q running inside the loop [6,26]: where N q c = 3, N l c = 1, Q f is the charge of the fermion, and B 1 is from eq. ( 6).
Decay widths into lepton pairs ℓ + ℓ − = e + e − , µ + µ − , τ + τ − are described by the formula The total width into leptons and the width into photons are shown in Fig. 7.

B. Hadronic decays
Let us now discuss the hadronic decays of ALPs.For m a ≫ Λ QCD , it is adequate to describe these decays by decays into quarks and gluons [26]: where c eff G is given by eq. ( 5), and N c = 3 is the number of colors.To approximately account for hadronization, when considering the kinematics, we replace the quark's mass with the mass of the lightest meson containing the given quark [50,51].For instance, for decays into cc, X q is a D meson.
Because of the same reason as it was discussed in the context of the Drell-Yan production (Sec.III C), decays into gluons dominate over decays into light quarks uū, d d, ss.However, above the D-meson pair production threshold, the decay into cc contributes a sizable fraction of the ALP total width because of a large c mass.
For m a = O(1 GeV), perturbative QCD breaks down, and one should use ChPT.To describe the palette of various decays (e.g., η/η ′ → ππγ, η ′ → 4π, or η ′ → ηππ) in agreement with the experimental data, the minimal ChPT is supplemented by phenomenological Lagrangians of the interactions of the pseudoscalar mesons P with vector [52,53], scalar [54], and tensor mesons [53], with the operators and their couplings being fixed by theoretical arguments (such as the chiral symmetry or anomaly matching conditions) and to match the experimental data on interactions of P .These mesons may contribute to the matrix elements as intermediate states; one example is the mixing of neutral vector mesons ρ 0 , ω, ϕ with photons.The ChPT width should match the parton-level width at some mass m a ≃ 1 GeV.
Ref. [10] followed this data-driven approach to describe the decays of the ALPs coupled to gluons; however, their approach also applies to the ALPs coupled to fermions.Ref. [41] repeated the analysis of [10] with some modifications for the ALPs with a non-universal explicitly isospin-breaking coupling to quarks.
In our analysis, we incorporate the ChPT Lagrangian (following the references above) in the Mathematica notebook accompanying the paper and calculate the matrix elements and decay widths for various processes (see Appendices B, A for details).We include the decay channels a → η2π, η ′ 2π, 4π, 3π, γ2π, 3η, ω2π, and 2V , where V = ω, K * , ϕ.As a cross-check, we reproduce the results of the SM decay widths of the mesons η and η ′ in the limit when the ALP matches them, i.e., when θ m 0 a = 1, f = f π , and m a = m m 0 .We have also qualitatively reproduced the results from [10] for the model of the ALPs coupled to gluons (see a discussion in Appendix A).
The summary of the hadronic widths for the ALPs is shown in Fig. 7.The decays of low-mass ALPs m a ≲ 1 GeV are saturated by a → 3π, a → γππ.At higher masses, decays into ηππ, 4π, and 2V become the dominant channels.The ChPT width matches with the width in the perturbative regime at m match ≃ 2 GeV.

C. Discussion
The leptonic, photonic, and hadronic widths are compared in Fig. 7.In Fig. 8, we show the branching ratios of the ALP decays into various final states.From For the ALP mass ma ≲ 2 GeV, ChPT describes decays of ALPs, while at higher masses, they may be approximated by the perturbative QCD; the vertical dashed line shows the matching mass.Bottom panel: the total non-hadronic (including di-e,µ,τ , and γ processes) and hadronic widths.The solid lines correspond to the scale Λ = 1 TeV, while the dashed ones correspond to Λ = mt.The vertical gray bands show the vicinity of ma = m π 0 /mη/m η ′ where the description via the ALP-meson mixing (and hence the results for the hadronic widths) breaks down.
the figures, we see that photonic decays are always subdominant, while hadronic widths are irrelevant for the phenomenology of light ALPs with m a ≲ 1 GeV, where leptonic decays dominate.For heavier ALPs, however, decays into hadrons dominate, increasing the total width by up to a factor of 100.This conclusion is in qualitative agreement with the paper [55], which studied a somewhat different model of a CP-odd scalar.We emphasize that in the mass range around 1 GeV, the hadronic decay width is significantly larger than the one obtained in Refs.[25,27] using perturbative QCD.The decay palette above 1 GeV is also qualitatively similar to the case of the ALPs coupled to gluons.This is because, for both of these models, the ALPs have mixing with the three neutral pseudoscalar mesons π 0 /η/η ′ .
Interestingly, the choice of the scale Λ practically does not influence the decay phenomenology.This is because it affects only the decays where the mixing with pions dominates among the others.These are the decays into 3π and γππ, which are important only in the mass range m a ≲ 1 GeV where leptonic decay widths are much larger.
To further stress the importance of the hadronic decays and finalize the discussion, we show the mass dependence of the ALP lifetime as computed in this work and the one widely considered in the past [8], when only leptonic The red and blue lines show, correspondingly, the lifetime in the approximation of using only leptonic width as in [8] and our result (both assuming Λ = 1 TeV), while the green line the results from [27] evaluated for the scale Λ = 4π TeV.
decays have been included.In the mass range 2m µ < m a < 2m τ , the full width in the description from [8] is saturated by dimuon decay.With hadronic decays being included, the branching ratio Br(a → µµ) is < 10% for masses m a ≳ 1 GeV.In Fig. 8, we show the branching ratios of the ALP decays into various final states (left panel) and the ALP lifetime (right panel).For the lifetime, we compare the predictions assuming the revised phenomenology and the description from past works: Ref. [8], which neglects hadronic decay modes and is widely used by the experiments community to derive constraints and sensitivities, and Ref. [27], which, following [25,26], approximates the hadronic decays by the decay a → 3π at m a ≲ 1 GeV and by the decays into a gluon and quark pairs above this mass.The lifetime from [8] is always much larger for m a ≳ 1 GeV.The lifetime from [27] coincides with our result for the ranges m a ≲ 1 GeV, m match < m a < 2m c , and m a ≳ 2m D .The origin of the discrepancy in the range 2m c < m a ≲ 4 GeV is that Ref. [27] turns on the decay a → cc above 2m c , even though this decay is kinematically impossible until the 2D threshold.
The previously neglected production and decay modes are expected to significantly change the landscape of the past constraints and future searches for ALPs.For example, let us consider searches for B → K(a →)µµ performed at LHCb [56,57].It is sensitive to the total ALP decay width as well as to the branching ratio Br(a → µµ).The constraints to ALPs are shown in Fig. 9, where, for comparison, we display the bound obtained assuming the ALP phenomenology description from Ref. [8] and the one obtained in this work.The updated constraints are weaker in the domain of large masses m a ≳ 1 GeV.Interestingly, for the revised phenomenology, the lower bound of the constraint lies in the regime where ALPs are shortlived and mainly decay within the detector, whereas for the old phenomenology, it mainly belongs to the parameter space of long-lived ALPs.We will revise further  [56,57] for the model of ALPs with the fermion coupling (3), assuming Λ = 1 TeV, for the plane ma − 2v h /f , where v h = 246 GeV is the Higgs VEV.The blue solid line: constraints assuming the ALP phenomenology obtained in this work.The red solid line: if assuming the phenomenology from [8], which only includes leptonic decays.The light lines of the same colors show the ALP decay lengths cτa⟨γa⟩ = 1 cm and 1 m assuming the corresponding phenomenology.
existing constraints, consider the ones previously not accounted for in the literature, and derive sensitivities for future experiments in a forthcoming work [58].

V. CONCLUSIONS
The model of ALPs universally coupled to fermions is considered by the physics beyond colliders (PBC) group as one of the benchmark models to test the potential of various experiments to explore the parameter space of feebly-interacting particles.Therefore, understanding the phenomenology of such ALPs is an important and timely question.In this work, we have revised the production and decay modes of ALPs at hadronic accelerator experiments, also considering the impact of the renormalization group (RG) flow of their couplings depending on the scale Λ at which the model is defined (see Sec. II and in particular Fig. 1).
For the production (Sec.III), we have considered decays of kaons and B mesons, the mixing with neutral mesons π 0 /η/η ′ , and the Drell-Yan process.For the production via mixing, we have found that the RG flow is very important, sizably changing the mixing angle squared between the ALPs and π 0 (Fig. 3).For the production from B mesons, we have included the decays B → X s + a, with X s being heavy kaon resonances K * 0 , K 1 , K 2 , . . ., which have not been considered previously in this context in the literature and increase the total production branching ratio by a factor of 3-4 for light ALPs (Fig. 2).Our results apply also to generic ALPs, provided that the low-energy Lagrangian describing the decay has the same operator expression.For the Drell-Yan production, we have considered the leadingorder fusion processes and shown that the cross-section suffers from a large systematic uncertainty (Fig. 4).
Depending on the scale Λ at which the ALP model ( 3) is defined and the collision energy, we have found that any of these processes may dominate the production (Fig. 5).
In particular, at DUNE collision energy, the production via mixing is the main production channel of the ALPs with mass below 2-3 GeV, while at larger masses, decays of B mesons are the main channel.At SPS energies, the hierarchy of production channels depends heavily on the scale Λ. Namely, if Λ is close to the EW scale, the main channels are mixing with mesons and the Drell-Yan process.Once Λ departs from the EW scale, decays via B mesons dominate.The main production channel may also depend on the geometric placement of the experiment (see Fig. 6).Finally, at the LHC, Λ practically does not influence the hierarchy and affects only the magnitude of the ALP flux.
We have studied the full palette of the ALP decays (Sec.IV), including the hadronic ones that were missing previously.In particular, we have found (Fig. 7 and also Fig. 8) that leptonic decays are the main channels for light ALPs with m a ≲ 1 GeV, while the hadronic decays dominate at higher masses, increasing the total ALP width by up to a factor of 100.Contrary to the production case, the decay widths are only weakly sensitive to the choice of Λ.
To simplify the use of our results by the community, we have implemented the ALP phenomenology studied in this work in a Mathematica notebook accompanying the paper.We have also implemented the model in SensCalc [35] -a public code evaluating sensitivities of different experiments.

FIG. 1 :
FIG.1:The RG flow of the quark and lepton couplings c ℓ , cq from the Lagrangian (3) (the left panel) and modulus of the flavor-changing abs coupling from eq. (9) at f = 1 PeV (the right panel).The scale Λ, at which the interactions (3) are defined (and where, by construction, c f (Λ) ≡ 1), is assumed to be between mt ≈ 173 GeV and 10 9 GeV.The values of the couplings are shown for the scale µ = 2 GeV.The vertical dashed line shows the reference scale Λ = 1 TeV used as our benchmark choice.

FIG. 3 :
FIG.3:The ALP mass dependence of the square of the mixing angle of the ALP with neutral mesons π 0 /η/η ′ (eq.(13)).The solid lines show the results assuming the scale Λ = 1 TeV, while the dashed line corresponds to Λ = mt.The vertical gray bands correspond to the vicinity of ma = m m 0 where the approximate description of the ALP interactions via the small mixing with the neutral mesons breaks down (see text for details).
FIG. 4: The probability of producing ALPs (3) in the DIS process (15) at various facilities.The bands show systematics uncertainties, which we obtained by varying the factorization and renormalization scales µ 2 r = µ 2 f = m 2 a by a factor of two (see text for details).

FIG. 5 :
FIG.5:The ALP production probabilities (eq.(16)) per proton collision at the LHC, SPS, and FNAL (DUNE), with the collision energies of √ s = 13 TeV, ≈ 28 GeV and ≈ 16 GeV correspondingly, for the model(3).The probabilities are evaluated for the value of the ALP coupling f = 1 PeV and assuming two scales Λ: Λ = mt (dashed lines) and Λ = 1 TeV (solid).The meaning of the vertical gray bands is the same as in Fig.3.See text for details.

FIG. 6 :
FIG.6:The production probabilities of the ALPs flying in the polar range θ < 10 mrad (solid lines) and θ > 10 mrad (dashed lines) at SPS, assuming the scale Λ = 300 GeV.By the curves "mixing", we summarize the production of the ALPs via the mixing angles with the mesons π 0 , η, and η ′ .The value f = 1 PeV is assumed.

f = 1 .FIG. 7 :
FIG.7: Summary of decays of the ALPs with the universal coupling to fermions(3).Top panel: hadronic decay widths, assuming the scale Λ = 1 TeV.For the ALP mass ma ≲ 2 GeV, ChPT describes decays of ALPs, while at higher masses, they may be approximated by the perturbative QCD; the vertical dashed line shows the matching mass.Bottom panel: the total non-hadronic (including di-e,µ,τ , and γ processes) and hadronic widths.The solid lines correspond to the scale Λ = 1 TeV, while the dashed ones correspond to Λ = mt.The vertical gray bands show the vicinity of ma = m π 0 /mη/m η ′ where the description via the ALP-meson mixing (and hence the results for the hadronic widths) breaks down.

FIG. 8 :
FIG.8: Left panel: branching ratios of the ALP decays into important final states (see Fig.7for the description of the states) for the model(3).Right panel: the ALP lifetime assuming f = 1 PeV.The red and blue lines show, correspondingly, the lifetime in the approximation of using only leptonic width as in[8] and our result (both assuming Λ = 1 TeV), while the green line the results from[27] evaluated for the scale Λ = 4π TeV.

FIG. 9 :
FIG.9: Re-interpretation of the model-independent LHCb constraints from the searches B → K + (FIP →)µµ reported in[56,57] for the model of ALPs with the fermion coupling (3), assuming Λ = 1 TeV, for the plane ma − 2v h /f , where v h = 246 GeV is the Higgs VEV.The blue solid line: constraints assuming the ALP phenomenology obtained in this work.The red solid line: if assuming the phenomenology from[8], which only includes leptonic decays.The light lines of the same colors show the ALP decay lengths cτa⟨γa⟩ = 1 cm and 1 m assuming the corresponding phenomenology.

TABLE I :
The values of the FCNC couplings from the Lagrangian (