Axion-Like Particles in Radiative Quarkonia Decays

Radiative quarkonia decays offer an ideal setting for probing Axion-Like Particle (ALP) interactions. This paper provides a comprehensive review of ALP production mechanisms through the $e^+ e^- \to \gamma\,a$ process at B- and Charm-factories, alongside an analysis of potential ALP decay channels. We derive constraints on ALP couplings to Standard Model (SM) fields, based on recent experimental results on quarkonia decays by the Belle II and BESIII collaborations. The analysis distinguishes between"invisible"and"visible"ALP decay scenarios. The"invisible"scenario, characterised by a mono-$\gamma$ plus missing-energy signature, enables stringent limits on ALP-photon and ALP-quark ($b$ or $c$) couplings. Moreover, extensive research at flavour factories has explored various"visible"ALP decays into SM final states, which depend on a larger set of ALP-SM couplings. To streamline the"visible"ALP scenario, we introduce additional theoretical assumptions, such as universal ALP-fermion couplings, or we adopt specific benchmark ALP models, aiming to minimise the number of independent variables in our analysis.


Introduction
Axion-Like Particles (ALPs) are a common feature of many extensions of the Standard Model (SM) of particle physics.Their lightness can be naturally justified if they are endowed with a pseudo-shift symmetry associated with the spontaneous breaking of an underlying global symmetry.A paradigmatic example is given by the QCD axion [1][2][3][4], which arises as a pseudo Nambu-Goldstone boson of a global U(1) PQ symmetry, anomalous under QCD and spontaneously broken at the scale f a ≫ v, where v denotes the electroweak scale.The main difference between the QCD axion and an ALP lies in abandoning the requirement that the only explicit breaking of the U(1) PQ symmetry arises from non-perturbative QCD effects, leading to the well known relation m a f a ≈ m π f π .Therefore, allowing the ALP mass, m a , and symmetry breaking scale, f a , to be independent parameters gives rise to a more general setup which is often described in terms of an effective Lagrangian containing operators up to d = 5 [5].The opportunity to look for ALPs with masses well above the MeV scale, whose couplings are not tightly constrained by astrophysical limits, opens up the possibility of probing ALP interactions at colliders [6][7][8][9][10][11][12][13][14][15] and via a broad class of rare processes, including flavour-violating [16][17][18][19][20][21][22][23][24][25][26][27] and CP-violating [28][29][30][31] observables.
Another relevant class of processes to probe ALP interactions is provided by radiative quarkonia decays of the type V → γa [32], with V = Υ, J/ψ a vector boson composed of heavy quarks.Historically, those observables played a fundamental role in ruling out the original Weinberg-Wilczek axion, see for instance Ref. [33].ALP searches at B-and Charm-factories span mass regions from a few MeV up to roughly 10 GeV.This ALP mass range is poorly constrained by astrophysical and cosmological probes, as well as by beam-dump experiments, making ALP searches via quarkonia decays essential in constraining the ALP parameter space.Moreover, recent new experimental results on quarkonia decays by the BESIII [34] and Belle II [35] collaborations motivate an updated and thorough phenomenological analysis, to constrain ALP couplings to the b and c quarks, alongside other ALP couplings, including those to photons.As argued e.g. in Ref. [36], the simultaneous presence of ALP couplings to quarks and photons, expected in general ALP setups, gives rise to new interesting phenomenological features.
In this work, we revisit the production of ALPs at B-and Charm-factories via the process e + e − → γa.Firstly, we update the Υ → γa analysis of Ref. [36], in which the ALP predominantly decays into an invisible channel (hereafter denoted as "invisible" ALP scenario), thus providing a mono-γ plus missing-energy signature, and extend the previous analysis to the charmonium sector.Moreover, we consider the complementary case in which the ALP decays visibly into the detector via SM final states.This scenario, denoted as "visible" ALP, opens more experimental channels to look for, albeit at the cost of enlarging the ALP parameter space.To reduce the number of independent parameters in the "visible" ALP scenario we make further theoretical assumptions, such as universal ALP-fermion couplings or consider a couple of benchmark ALP models, inspired by standard QCD axion models.
The paper is structured as follows.In Sec. 2 we discuss the general form of the ALP effective Lagrangian below the electroweak scale and introduce two benchmark ALP models that will be used in the following phenomenological analysis.In Sec. 3 we provide the general framework for ALP production in the context of B-and Charm-factories, while in Sec. 4 we collect general formulae for ALP decays into SM final states.Sec. 5 is devoted instead to the phenomenological analysis, setting constraints both on the "invisible" and "visible" ALP scenarios.Our results are summarised in Sec.6, while Appendix A focuses on a technical issue, that is the exact diagonalisation of the ALP-pion mixing, which is relevant for our analysis and could also hold broader interest.

The ALP Lagrangian
The most general CP-conserving and flavour-diagonal effective Lagrangian, describing the ALP-SM particle interactions at energies below the electroweak scale is given by the following d = 5 operators: with f running over all the SM fermions but the top-quark and V µν ≡ ϵ µναβ V αβ /2 (with ϵ 0123 = 1) the dual field strength.Note that when considering ALP scenarios one typically works in the limit in which the ALP decay constant, f a , is much larger than the ALP mass, i.e. f a ≫ m a .All couplings in Eq. (2.1) are defined at the scale of the experiment, i.e. c aXX ≡ c aXX (µ = m Υ , m J/ψ ).We here assume that the heavy SM fields, i.e. top, Higgs and weak gauge bosons, have been integrated out and their effects are already encoded in the low-energy couplings by matching the high-and low-energy effective Lagrangians at the electroweak scale, and then running the Wilson coefficients down to the scale of the experiment (for details, see Ref. [37]).Another commonly used notation, especially in experimental papers, incorporates the scale f a in the coupling's definition and introduces the following mass-dimensional couplings: In particular, notice the different α em (α s ) normalization of the ALP-photon (gluon) coupling in the two versions of the Lagrangian, with the loop-origin of the photon (gluon) coupling being highlighted in the notation of Eq. (2.1).ALPs may also act as portals to a light-dark sector.In this case, additional ALP couplings are introduced.Assuming, for simplicity, the existence of a single light and neutral dark fermion χ, the following term is added to the d = 5 ALP effective Lagrangian: with c aχχ a coupling that can induce a sizeable ALP decay into invisible final states, whenever m a ≳ 2 m χ .

Benchmark ALP models
The effective field theory approach of Eq. (2.1) features several Wilson coefficients, which do not allow in general for a simple representation of the experimental constraints from quarkonia decays.Therefore, it is also useful to consider ultraviolet (UV) complete scenarios which provide nontrivial correlations among the Wilson coefficient, thus considerably reducing the ALP parameter space.In the following, we will provide a couple of benchmark ALP frameworks, inspired by the canonical DFSZ [38,39] and KSVZ [40,41] axion models.Differently from the case of the QCD axion, we will assume here that the ALP mass is a free parameter.See also Ref. [42] for related examples.

DFSZ-QED ALP
A useful benchmark scenario is provided by a variant of the DFSZ model with no QCD-PQ anomalous couplings, that we denote as DFSZ-QED.This is obtained by means of the following Yukawa Lagrangian, with the two Higgs doublets denoted as H q and H ℓ The scalar potential is assumed to contain an operator of the type H † q H ℓ ϕ †2 , with ϕ a SMsinglet complex scalar.This operator constrains the PQ charges of the scalar fields, generically denoted by X q,ℓ,ϕ and normalised so that X ϕ = 1, to satisfy X ℓ − X q = 2X ϕ = 2. Imposing the orthogonality of the PQ and hypercharge currents and defining tan β ≡ t β = v ℓ /v q , with v 2 ew ≡ v 2 q + v 2 ℓ = (246 GeV) 2 , one obtains X q = −2s 2 β and X ℓ = 2c 2 β .It can be readily verified that this model implies no QCD-PQ anomaly since SM quarks interact with the same Higgs doublet.The ALP couplings to SM fields can then be obtained following standard techniques (see e.g.[43]) and, matching onto the notation of Eq. (2.1), one finds DFSZ-QED: where the ALP-SM fermions couplings are understood to be universal.

KSVZ ALP
As a second benchmark setup, we consider a KSVZ-like ALP model in which the PQ anomaly is carried by a new coloured heavy fermion, Q, charged under the PQ symmetry, with Yukawa interaction Q L Q R ϕ, and ϕ a SM singlet field carrying unit PQ charge.After integrating out the new heavy fermion, assuming e.g. that it transforms in the fundamental of colour, one obtains (see e.g.[43]) for the ALP-SM tree-level couplings.

ALP production mechanisms
The ALP production mechanisms at B-and Charm-factories via the process e + e − → γa have been carefully analysed in Ref. [36], to which we refer for most of the details.In the following, we will shortly recall the main features that will be relevant for the phenomenological discussion of the subsequent sections.The simplest way of producing an ALP in e + e − colliders is via the non-resonant tree-level process e + e − → γa, that is dominated by the s-channel photon mediated diagram.The total non-resonant cross section, in the centre-of-mass rest frame, is given by [28]: As we are interested in a relatively large ALP mass range (from few MeV up to about 10 GeV) the m a dependence is fully accounted for in Eq. (3.1), while terms proportional to m2 e /s can be safely neglected at the energies relevant for flavour factories.Note that the non-resonant ALP-photon production via a t-channel electron contribution is usually neglected as suppressed by the electron mass: the underlying assumption here is that c aee /c aγγ ≲ 10 2 (or equivalently g aee /g aγγ ≲ 10 4 ). 1hile the non-resonant contribution to ALP production in Eq. (3.1) is unavoidable in any e + e − experiment, the situation at B-and Charm-factories is more intricate since these experiments operate around specific resonances (i.e.Υ(nS) or J/ψ(nS)).It is, therefore, crucial to correctly account for the resonantly enhanced contributions.Generalising the discussion in [36], the resonant cross section for an ALP production via vector quarkonia resonances, V = Υ, J/ψ, in the Breit-Wigner approximation reads: where m V and Γ V are, respectively, the mass and the width of the specific resonance and σ peak is the peak cross section defined in terms of the leptonic branching fraction, B(V → e + e − ), experimentally determined for each different intermediate state [44].The effective ALP couplings defined in Eq. (2.1) enter in the B(V → aγ) branching fraction as [36]: with Q Q and c aQQ the electromagnetic charge and the ALP-fermion couplings of the quarkonia valence quarks, i.e.Q = c, b for the cases under consideration 2 , while f V are the quarkonium decay constants, which can be obtained from [45].Naively, one would expect the resonant contribution to always dominate by orders of magnitude over the non-resonant one.However, this expectation is incorrect in some of the cases under consideration, since the resonance widths are typically much narrower than the energy beam resolution, σ W .Therefore, when the resonance is not fully resolved by the experiment, a convolution of the theoretical resonant contribution with a Gaussian spread function has to be performed (see e.g.[46]) and an "experimental" resonant cross section can be defined as: In the case of very narrow resonances, the previous result can be simplified to: with the parameter ρ accounting for the suppression of the "experimental" resonant cross section due to the non-negligible beam-energy spread.
Depending on the experimental setup, ALP searches at flavour factories can be classified into three different categories, based on the different ALP production mechanisms.i) Non-resonant searches: This scenario arises when the ALP is produced off-resonance or the resonance is very spread, as for example is the case of Belle II searches at the Υ(4S) resonance [35].These searches are dubbed non-resonant as in this case the non-resonant cross section in Eq. (3.1) has to be used and therefore only a sensitivity to g aγγ can be obtained from the ALP production.
ii) Resonant searches: Excited quarkonia states can decay into lighter quarkonia resonances via pion emission, for example, Υ(2S) → Υ(1S)π + π − .By exploiting the kinematics of the final states, the lighter meson resonance can be fully reconstructed and its mono-γ decay mode Υ(1S) → γ a analysed, as done for example in Ref. [47].These searches are dubbed in the following as resonant since they allow to directly probe B(V → γa) in Eq. (3.4) and give information on a specific combination of g aγγ and g aQQ couplings.
iii) Mixed searches: Experimental searches sometimes are performed at an energy that is around a specific (narrow) resonance √ s ≈ m V (nS) , but without identifying it kinematically.
From the above discussion it is clear that in this case, due to the larger beam energy uncertainties, the production mechanism depends on both the resonant and non-resonant contributions, being sensitive to a different g aγγ and g aQQ combination, compared to the resonant searches.In the case under examination, the following simple prescription provides a good approximation of the production cross section [36]: In Tab. 1 an estimate of the "experimental" resonant cross sections compared to the non-resonant ones is presented, for the relevant decay channels.The experimental energy spread of the e + e − beams at flavour factories can be typically taken in the range σ W ≈ 2 − 5 MeV [44,[48][49][50].This value is considerably larger than all the considered resonance widths, typically in the 20 − 100 keV range, with the only exception of the Υ(4S) resonance, for which Γ Υ (4S) = 20.5 MeV.As it can be seen, for the Υ(4S) case the Gaussian smearing is practically ineffective, i.e. ρ ≈ 1 being Γ Υ (4S) ≫ σ W .However, at the same time, the broadness of the resonance largely reduces the resonant contribution, five orders of magnitude below the non-resonant one, due to the smallness of the corresponding σ peak value.In all the other cases, instead, the smearing procedure strongly suppresses the naive Breit-Wigner theoretical expectation by roughly a factor of 10 −3 .In these cases, the "experimental" resonant cross section is smaller than the non-resonant one, even if it can still contribute with numerically significant effects between 20% and 50% of the non-resonant one, which should therefore be taken into account when interpreting experimental searches.For illustrative purposes, the numbers reported in Tab. 1 have been calculated assuming vanishing ALP-fermion couplings and, as a consequence, the c aγγ dependence simply cancels out in the cross section ratio.In the general case, with more ALP couplings, results can slightly change.
A detailed discussion regarding all the subtleties entering in the interpretation of the e + e − → γ a production at the Υ(nS) resonances has been performed in Ref. [36], to which the interested reader is referred.

ALP decay channels
Being interested in ALP production from radiative quarkonia decays, in the following, we will focus on ALP masses m a ≲ 2 m q , with q = c, b depending on whether V = J/Ψ or Υ.The total ALP decay width in SM particles reads: where Θ(x) is the Heaviside step function, N f c = 1(3) for leptons (quarks) and f = {e, µ, τ, c, b} extends to all kinematically accessible fermions, except for light quark flavours which are included in the the light hadronic decay width, Γ a→lh , that is provided in Eq. (4.2), while the effective ALP couplings to photons is defined below in Eq. (4.6).
In the m a ∼ O(GeV) region, ALP decays into light hadrons become at the same time relevant and difficult to compute, requiring different approaches depending on the value of m a .While for m a ≲ 1 GeV one can use chiral perturbation theory to predict exclusive hadronic ALP decays, for m a ≳ 2 GeV one can rely on the quark-hadron duality [51,52] to compute the inclusive ALP decay rate into hadrons.However, in the m a ∈ [1, 2] GeV region, one lays outside the applicability range of both chiral perturbation theory and perturbative QCD.For this reason, the 1-2 GeV ALP mass range will be conservatively excluded when deriving bounds on the ALP-SM couplings 3 .The ALP decay width into light hadrons, in the two m a regions considered, is then given by [9,54]: where n q is the number of active quark flavours and we have introduced the ALP-pion and ALP-gluon effective couplings, defined respectively by with the loop function B 1 (τ ) given by respectively below and above the quark production threshold.The terms proportional to the functions g 00 and g +− in Eq. ( 4.2) correspond respectively to the three-body decays of a → 3π 0 and a → π + π − π 0 , while other three-body decays containing electrons or photons are suppressed by powers of α em .The explicit expression of the functions g 00 and g +− can be found in Ref. [9], where the reader is referred for additional details.It should be mentioned that the chiral description of a → 3π processes breaks down just above kinematical threshold [55], so although we expect these decay widths to be O(1) correct, a more refined description of the ALP width below the GeV scale would require techniques beyond standard chiral perturbation theory (see e.g.[56] for a related example).Finally, in Eq. (4.1) the effective photon coupling, c eff aγγ , has been defined to include all possible one-loop contributions.Again, this coupling needs to be defined differently in the two considered ALP mass regions.Neglecting the subleading top-quark and W loops one obtains [9,37] with Q f the fermions' electromagnetic charge.An approximate expression for the ALP-pion mixing term, U aπ , valid for More details on the ALP-pion diagonalisation procedure and the exact expression for the mixing angle are reported for completeness in Appendix A.
In Fig. 1 we display the branching ratio B(a → SM), based on Eq. (4.1) and assuming, for the sake of illustration, all the Wilson coefficients to be c aXX = 1 and f a = 1 TeV.In this way, one  can compare the relative strengths of the ALP decay to SM final states.In the sub-GeV ALP mass region, the ALP decays preferably in electrons and muons, except in the proximity of the maximal mixing region, i.e. m a = m π , where the ALP-photon coupling gets largely enhanced, while for m a ≳ 2 GeV the larger ALP branching ratio is dominantly into light and heavy hadrons with the τ contribution is only slightly subdominant when kinematically allowed.Note, however, that in the conventions of Eq. (2.1), the fine-structure constant α em is factorised outside the ALP-photon coupling, c aγγ .Nonetheless, the non-decoupling nature of the anomalous terms may easily give rise to larger contributions, depending on the multiplicity of heavy electromagneticallycharged exotic fermions in the UV-completed model (see e.g.[57,58]).If, instead, the alternative definitions of Eq. (2.2) for the ALP-SM couplings are used, and g aXX = 1 is assumed, the ALPphoton branching ratio would be enhanced by a factor ∼ 10 5 , thus resulting in the dominant contribution in most of the considered m a range.Finally, in a more general setup, one can consider the possibility in which the ALP is the portal to a light-dark sector.In the following, we will consider the simplest scenario described by the Lagrangian of Eq. (2.3), where a single light exotic fermion has been introduced with coupling c aχχ .Then, in all generality, the cross section for a specific ALP decay into the final state where σ prod (s) is the relevant production cross section of the quarkonia state (resonant, nonresonant or mixed) and B(a → X) the branching fraction defined as Two simplified scenarios for the ALP decay are typically envisaged in the literature, dubbed respectively "invisible" and "visible".For the "visible" ALP scenario we assume c aχχ ≈ 0 and therefore B(a → SM) ≈ 1.Then the relative strength of all SM contributions is the one already Exp.
discussed and depicted in Fig. 1.Conversely, the "invisible" ALP scenario is obtained by imposing a large hierarchy between the ALP couplings to the SM and dark sector particles, i.e. c aχχ ≫ c aSM , which therefore leads to B(a → SM) ≈ 0 and B(a → χχ) ≈ 1, within our simplified framework.
5 Phenomenology of the "invisible" and "visible" ALP scenarios In this section, we examine the phenomenological analysis of the J/Ψ and Υ radiative quarkonia decays, addressing the "invisible" and "visible" ALP scenarios separately.

The "invisible" ALP scenario
In this case, the ALP decays invisibly inside the detector, with a very clean missing-energy signature and typically well-controlled SM backgrounds.In the specific case of radiative quarkonia decays, the main experimental signature is simply an isolated photon plus missing energy.This scenario is the most economical in terms of fitting parameters as it depends only on the ALP couplings g aγγ and g aQQ (with Q = c or b) through the production cross sections in Eqs.(3.1) and (3.4).A detailed analysis of the invisible ALP decay at the Υ(nS) resonances can be found in Ref. [36], while a more general study of hadronic and leptonic meson decays in invisible ALPs can be found in Refs.[22,23,26].A summary of the experimental searches considered in the "invisible" ALP analysis is provided in Tab. 2. Although there are no new experimental data for the Υ(nS) radiative decays into an "invisible" ALP, with respect to the analysis provided in [36], BESIII has performed in [62] a search for a J/Ψ radiative decay through the process Ψ(3686) → π + π − J/Ψ, followed by the decay J/Ψ → γ + invisible, assuming in the final state an exotic neutral scalar decaying invisibly inside the detector.As the J/Ψ is kinematically reconstructed from the original Ψ(3686) decay process, one can assume the ALP being produced with the fully resonant cross section of Eq. (3.2), thus acquiring a sensitivity to the (c aγγ , c acc ) couplings for an ALP mass up to m a ≲ 1.2 GeV.This analysis can complement B-factories limits on the ALP-photon coupling in the lower part of the ALP mass range, and provides new direct bounds on the ALP-charm interaction.In Fig. 2, these new bounds are shown in red and compared with previous limits obtained in [36], depicted as blue and green dashed lines.
In addition to the invisible searches mentioned above, valuable information can also be obtained through the recast of recently published limits on the quarkonia ALP decay rates in the di-photon channel, e + e − → V → γ (a → γγ) from BESIII [34] and Belle II [35] collaborations, highlighting expectations from future "invisible" ALP searches at both these facilities.The BE-  SIII search of Ref. [34] looks for an ALP through the process Ψ(3686) → π + π − J/Ψ, followed by the resonant decay J/Ψ → γ (a → γγ).In the BESIII analysis only the ALP-photon coupling is considered and consequently the ALP-SM branching ratio B(a → γγ) = 1 is assumed.Using Eqs.(3.2)-(3.4),and setting c acc = 0, from the BESIII limits on B(V → γa) one can derive bounds on the ALP-photon coupling.In the left panel of Fig. 2 our exclusion limits on g aγγ are shown as light-red region, perfectly in agreement with the bounds presented in [34].
Notice, however, that BESIII limits on B(V → γa) obtained from the tri-gamma channel can also be used to provide bounds on the foreseen sensitivities of the collaboration on the ALP-SM couplings in the "invisible" ALP scenario, to be compared with the (g aγγ , g aQQ ) limits derived from BESIII invisible, Belle and BaBar searches.From the two plots in Fig. 2 one can see that a future BESIII invisible ALP search would feature a sensitivity, up to the kinematically available mass region, fully comparable with B-factories limits.It should be stressed that such a limit on the g acc coupling, once provided, would represent the most stringent one in the considered ALP mass region, with B-factories yielding a comparable bound on the ALP-bottom coupling.
On the same footing, Belle II has recently published in Ref. [35] a similar analysis looking for an ALP through the process e + e − → Υ(4S) → γ (a → γγ), thus providing bounds on the ALP-photon coupling, once B(a → γγ) = 1 is assumed.As explained before, the Υ(4S) resonance is so spread that the non-resonant ALP production cross section largely dominates and Eq.(3.1) has to be used in order to derive bounds in the ALP-photon coupling.Again, as before, the obtained bounds on the ALP non-resonant cross section can also be recast as bounds on the g aγγ couplings in the "invisible" ALP scenario.The resulting exclusion limits are shown as a dark grey region in the left panel of Fig. 2, while clearly in this search no sensitivity can be obtained on the ALP-quark coupling.
Finally, we present in Fig. 3 all the available information on the "invisible" ALP scenario from J/Ψ and Υ radiative decays in the (g aγγ , g aQQ ) parameter space, with Q either c or b depending on g aQQ /g aγγ < 0 m a = 0.2 GeV m a = 2 GeV Figure 3: Excluded (g aγγ , g aQQ ) parameter space in the "invisible" scenario, for two different values of the ALP mass m a = 0.2 GeV (left) and m a = 2 GeV (right) and for the two possible choices of the sign of the coupling ratio, g aQQ /g aγγ > 0 (upper) and g aQQ /g aγγ < 0 (lower).Q = c or b respectively for Charm-or B-factories.Dashed constraints are taken from [36].The dark-red area is from the BESIII search [34], while the continuous black line displays the Belle II search [35].
the corresponding experiment and for two different reference values of the ALP mass, m a = 0.2 (left) and 2 GeV(right).The dark and light red areas show the (g aγγ , g acc ) exclusion regions from the invisible and tri-gamma (recast) BESIII analyses, while the blue and green dashed line indicate the bounds on (g aγγ , g abb ) from BaBar and Belle Υ Invisible radiative decays (taken from Ref. [36]).The full black line shows, instead, the bound obtained from the Belle II tri-gamma (recast) analysis.In the two upper plots, one can observe the flat direction associated with resonant searches when g aQQ /g aγγ > 0 is chosen, absent instead when the opposite sign choice is performed (lower plots).This otherwise unconstrained direction in the (g aγγ , g aQQ ) parameter space can be resolved only by means of the BaBar Υ(3S) mixed search (blue dashed line) or the Exp.
The searches are classified as resonant, non-resonant and mixed, according to the experimental setup.
Belle II non-resonant one (black continuous line).
As a final comment, notice that, at least in the sub-GeV mass range, BESIII has already the same potential sensitivity both on the ALP-photon and ALP-quark couplings than B-factories.This highlights the importance of undertaking a real "invisible" ALP search, updating the 2020 result of Ref. [62] by using the present luminosity and extending the search up to m a = 3 GeV.

The "visible" ALP scenario
Let us consider now the case B(a → SM) ≈ 1, so that the ALP decays visibly into SM final states inside the detector.This scenario, dubbed "visible" ALP, opens up more experimental channels to look for, albeit at the cost of enlarging the ALP parameter space.
Several experiments have looked at different decay channels in the allowed ALP mass region.The most significant searches are summarised in Tab. 3, where we have explicitly indicated also the ALP production mechanism, as described in Sec. 3, and the associated decay mode.The oldest data come from the BaBar experiment, where searches in several decay channels have been released.For this analysis, we make use of the data from the Υ(1S) radiative decays [64,65,68] into a scalar particle, that subsequently decays visibly into muons and c-quarks.In these searches the Υ(1S) is kinematically reconstructed from a higher quarkonium resonance decay and hence the contribution is resonant.In addition, we have also used the BaBar dataset of Υ(3S) and Υ(2S) decays into light and charm hadrons [63].Since the latter are untagged, the mixed ALP production cross section has been used in this case.Recall, however, that the ALP decay into light hadrons can only be approximately estimated (see discussion in Sec. 4) and hence the bounds derived through these data should be interpreted with some caution, especially in the sub-GeV ALP mass range.Previous Belle visible searches [66] with the ALP decaying into muons and taus are also considered in this work.The latest available analysis of radiative Υ decay into a visible ALP is the Belle II measurement of Ref. [35] of the Υ(4S) decay into 3γ.As explained in Sec. 3, in this case, the non-resonant ALP production cross section has to be used, due to the very spread nature of this specific resonance.Newer data coming from BESIII complement B-factories searches with the ALP production via the radiative J/Ψ decay and the subsequent ALP decay into muons [67] and photons [34].These data are important because they can provide additional direct information on the ALP-charm coupling, inaccessible at B-factories.
Due to the large number of independent parameters entering in the "visible" ALP decay scenario, in the following, we will present our phenomenological analysis in two simplified frameworks, considering first universal ALP-fermion couplings and then employing the two benchmark ALP models introduced in Sec.2.1.
Figure 4: BESIII, Belle(II) and BaBar 90% exclusion bounds obtained from visible ALP decay searches, in the universal ALP-fermion scenario and assuming vanishing tree-level ALP-gauge anomalous couplings.

Universal ALP-fermion coupling
As a first simplified scenario, we study how available experimental data from radiative quarkonia decays can constraint the ALP parameter space by assuming a universal ALP-fermion coupling, denoted hereafter as c af f .This universal ALP-fermion coupling scenario can be theoretically motivated by generating all the ALP-fermion couplings exclusively via the ALP-Higgs operator (∂ µ a)H † ↔ D µ H, see for instance Ref. [5].We neglect here, for simplicity, the small differences in the individual ALP-fermion couplings induced by running effects.
In Fig. 4 the bounds on the universal ALP-fermion coupling, g af f ≡ c af f /f a , obtained from radiative quarkonia decays into a visible ALP are shown, assuming vanishing tree-level ALPgauge anomalous Wilson coefficients, c agg = c aγγ = 0.However, for consistency we are including the fermion loop-induced contributions, proportional to c af f , as discussed in Sec. 4. In Fig. 4 one can see that for m a ≲ 1 GeV muon searches at Charm and B-factories give very similar results, g af f ≲ 3 × 10 −3 GeV −1 at 90% CL.No relevant constraints can be obtained in the ALP lowmass region from ALP hadronic decays as the ALP decay width into light-hadrons of Eq. (4.2) is almost vanishing due to an "accidental" cancellation that takes place in the ALP-fermion universal scenario, i.e. c auu = c add , leading to c aπ ≈ 0, see Eq. (4.3).Conversely, at higher ALP masses, m a ≳ 2 GeV, the strongest bounds on the universal ALP-fermion coupling come from BaBar searches in ALP-hadron (yellow line) and ALP-tau (green line) decays when kinematically allowed, with all the other channels providing bounds 5-10 times weaker and with BESIII data still subleading compared to the same BaBar/Belle channel.
Bounds on g aγγ for the "visible" ALP scenario, assuming a vanishing universal fermion cou- pling, g af f = 0, are not shown here as they coincide with the excluded light-pink and grey regions depicted in the left plot of Fig. 2, obtained from the BESIII and Belle II 3γ searches of Ref. [34,35].
Finally in Fig. 5 the combined bounds on (g aγγ , g af f ) from the "visible" ALP searches are summarised, assuming a vanishing tree-level ALP-gluon coupling, g agg = 0.In the left plots, the excluded parameter space is shown for a value of the ALP mass well inside the region of validity of chiral perturbation theory, i.e. m a = 0.25 GeV.In these plots, for the sake of simplicity, we choose to display the BESIII muons data alone, as all muons constraints practically overlap, as learned from Fig. 4. In the upper/lower plots the bounds are shown for the two alternative sign choices, g af f /g aγγ ≶ 0. In the upper-left plot, one can observe that the typical flat direction present in all resonant channels (red line) is resolved by the mixed BESIII muon channel and by the nonresonant 3γ Belle II search at the Υ(4S) resonance.Therefore, once again the complementarity of searches with a different ALP production mechanism appears evident.Clearly, no flat direction appears when the negative sign is chosen.
In the right plots of Fig. 5 similar bounds are shown but for an ALP mass well inside the perturbative QCD validity range, i.e. m a = 5 GeV.For this large ALP mass, only data from B-factories are available.The resonant data from tau and c-quark ALP decays give constraints slightly weaker than the mixed ALP-hadrons decay channel.However, all these experiments give relevant bounds only on the ALP-fermion coupling while the non-resonant Belle II 3γ search becomes fundamental in closing the two-dimensional parameter space.In this case, the tau and hadronic searches are the most relevant ones, setting a robust bound on the universal coupling g af f ≲ 1 − 2 × 10 −3 GeV −1 .Finally, it can be seen that bounds on the fermion coupling are an order of magnitude better in the low-mass scenario than in the heavier ALP case.This fact can already be noticed from looking at the individual coupling bound of Fig. 4, and similarly for the photon coupling in Fig. 2, and it is associated to a phase space suppression at higher ALP masses.

Benchmark ALP models
Another possibility, in order to reduce the number of free parameters in the "visible" ALP scenario, is to consider the UV-complete ALP models (previously introduce in Sec.2.1) which provide non-trivial correlations among the Wilson coefficient, thus considerably reducing the ALP parameter space.Notice that in the effective ALP Lagrangian of Eq. (2.1) we have assumed that top-quark (and the heavy weak bosons) contributions were already matched into the low-energy Wilson coefficients.However, when considering specific UV-complete models these contributions need to be included explicitly.To do so, we implement in our numerical analysis the analytical results of Ref. [24] and take into account the effects of running from the high-energy scale f a (1 − 10 TeV in our benchmark cases) and integrating out the top at the electroweak scale.This procedure can account up to a 10% modification of the UV couplings defined in Eq. (2.5) and (2.6).On the other hand, the weak boson contributions are safely negligible [24].
The first model we have considered is the DFSZ-QED ALP, which is anomalous under electromagnetism but not under QCD.It is not completely fermion universal as all ALP-quarks couplings are the same but different from the ALP-lepton ones.In the DFSZ-QED model, there are three independent parameters, m a , f a and tan β.In the following, we will set for simplicity tan β = 1, as in this point both couplings to leptons and quarks have the same magnitude.For this benchmark point, we plot in the left panel of Fig. 6 the excluded 1/f a region as a function of the ALP mass m a .Notice that photon searches are dominant only when the ALP decay into muons is not kinetically allowed.The BESIII bound is sizeable because the production of the ALP is made via the c coupling, while the 3γ Belle II search is of non-resonant type and hence is only sensitive to the photon coupling suppressed by α 2 em in our conventions of Eq. (2.1).When the ALP decay into muons is allowed this channel becomes the most stringent one below 1 GeV, with all experiments having similar sensitivities.At higher masses, hadronic and tau decays set the most stringent bounds, while cc exclusive decays are less relevant.As all fermion couplings are generated at tree level in this model the bounds are similar to those of the universal scenario shown in Fig. 4, ranging from f a ∼ 10 TeV for lower masses and f a ∼ 1 TeV in the heavier mass region.
In the second ALP model considered, the KSVZ ALP, all couplings are generated via the c agg anomalous term: including one-loop-induced quark couplings (see Ref. [42]), and a photon coupling via mixing with the pion.In this case, one has only two independent parameters (m a , f a ), and expects in general weaker constraints on the scale f a .We do not include in this study the two-loop ALP-lepton and ALP-photon contributions, as they are safely negligible.The obtained bounds on 1/f a as a function of the ALP mass are shown in the right panel of Fig. 6.For low values of m a the 3γ BESIII search is the most relevant, together with the 3π channel, once kinematically allowed.Limits from Belle II are very weak, as the production is only via the photon coupling and has therefore an extra α em suppression.Similarly to the DFSZ-QED ALP benchmark, at larger masses, the strongest bound is set by the hadronic decays, while the limits from exclusive c−quark searches are still an order of magnitude weaker.

Conclusions
In this paper we have revisited the ALP production in Υ and J/ψ radiative decays.The ALP production can proceed through three distinct and complementary experimental setups, dubbed as non-resonant, resonant and mixed.Depending on the specific production mechanism, they enable to access different combinations of ALP-SM particle couplings.We have then considered both the "invisible" and "visible" ALP decay scenarios, depending respectively on whether the ALP escapes detection (e.g. by decaying into an unspecified dark sector) or decays into SM final states.
For the "invisible" ALP scenario, we have updated the analysis presented in Ref. [36] by adding the constraints on the (g aγγ , g acc ) parameter space from the J/ψ → γ + E mis BESIII analysis [62].In addition, we presented here the projected sensitivity on the (g aγγ , g acc ) and (g aγγ , g abb ) parameter space of future "invisible" ALP searches, obtained by recasting the recent BESIII J/Ψ → γa (a → γγ) search [34] along with the constraints from the analogous Belle II Υ(4S) → γa (a → γγ) channel [35].ALP-photon and ALP-Q couplings, with Q = c, b, larger than roughly 10 −3 GeV −1 are already excluded by present data for ALP masses in the m a ∈ [0.1, 10] GeV range.Assuming the already achieved BESIII and Belle II luminosities, a factor of five improvement in the above limits is expected when the new analyses in the mono-γ plus missing-energy channel will become available.
In the case of a "visible" ALP decay, two m a regions have been separately considered, respectively m a ≲ 1 GeV and m a ≳ 2 GeV.This is because the ALP decay width, Γ(a → SM), turns out to be difficult to estimate in the intermediate 1-2 GeV region, where neither chiral perturbation theory nor perturbative QCD can be applied.The landscape of currently available searches for visible ALP decay channels from radiative quarkonia decays has been then presented.To streamline the analysis of the "visible" ALP scenario, we have introduced additional theoretical assumptions, such as "universal ALP-fermion couplings" or specific "benchmark ALP models", aimed to reduce the number of independent parameters involved.As a general feature, we noticed that the strongest constraints on the ALP-SM couplings in the sub-GeV m a range come from the BESIII and Belle II muon and 3γ channels, with an almost equivalent sensitivity on f a ≳ 10 3 − 10 4 GeV.For the multi-GeV m a case the bounds are typically one order of magnitude weaker than for the lower ALP mass range, with the strongest constraints coming from hadronic channels at B-factories.
From the analysis presented in this paper clearly emerges the relevance of quarkonia decays in constraining the ALP parameter space.Further significant improvements on the ALP-fermion flavour-conserving couplings can be envisaged when/if future experimental searches at B-and Charm-factories will be available.In particular, BESIII searches in the light hadrons channel could be extremely useful to further constrain the ALP coupling to pions.In the multi-GeV ALP range, instead, Belle II radiative decays of Υ(4S) into both invisible and visible channels, besides the already available 3γ channel, would offer a complementary set of information thanks to the underlying non-resonant ALP production mechanism.Finally, quarkonia decays are also very promising for constraining ALP-fermion flavour-violating couplings, which were not addressed in this study but will be the focus of a future work.
In order to implement the transformation that diagonalizes the mass matrix on the ALP and pion fields, it is actually more convenient to provide an expression directly for the trigonometric functions: Here, the signs should be properly chosen according to the discussion above (cf.also Fig. 7).
We then conclude that the outer signs in Eqs.(A.12)-(A.13)need to be both positive.On the other hand, for ϵ > 0 and m π < m a , we have 3π/4 < θ < π, and hence we need to choose plus for the sine in Eq. (A.12) and minus for the cosine in Eq. (A.13).Note that the sign of ϵ will determine whether the sine is larger or smaller than the cosine.It is also useful to note that it is equivalent to diagonalise doing a θ + π rotation.This is helpful if one wants to display a dependence of a certain observable with respect to m a , as if we diagonalise in the case ϵ > 0 and cross m a = m π then one needs to switch the sign of the cosine in a non-continuous way.
Then the full diagonalization matrix acting on the physical mass eigenstates can be defined as are the physical states associated to the eigenvalues m 1 and m 2 .We will identify these states with what we call as the "physical" ALP and pion fields later on, as this requires one last argument.
Away from the m a ≈ m π region, the transformation connecting the a and π to the physical fields, is usually approximated as which are valid for m 2 π − m 2 a /max{m 2 π , m 2 a } ≫ f π /f a .Otherwise, the full formulae in terms of U aπ should be employed.
As a final remark, note that if we take the limit f π /f a → 0 then Eq. (A.17This limit explains why the branching ratio into photons is independent of ϵ.This is reflected in Fig. 1, where modifying the value of f a does not change the shape of the branching ratio to photons as long as f π ≪ f a , making it independent of ϵ.

Figure 1 :
Figure 1: Branching ratios of ALP decays into SM particles as a function of the ALP mass, setting all Wilson coefficients c aXX = 1 and f a = 1 TeV.

Figure 2 :
Figure 2: Limits from the BESIII invisible (red) and tri-gamma recast (dark red) searches in the "invisible" ALP scenario.The dashed blue, green and black lines represent respectively BaBar, Belle II invisible and the Belle II tri-gamma recast limits and are shown for comparison.In the right panel, Q = c or b respectively for Charm-or B-factories.

Figure 5 :
Figure5: Excluded parameter space for the two benchmark values m a = 0.25 GeV (left) and m a = 5 GeV (right) and for the two possible choices of the sign of the coupling ratio, g af f /g aγγ > 0 (upper) and g af f /g aγγ < 0 (lower).The blue line indicates the BESIII search of an ALP decaying into two muons, which largely overlaps with the BaBar and Belle searches.

Figure 6 :
Figure 6: Bounds on 1/f a as function of the ALP mass, m a , for the two benchmark ALP models considered in the text: DFSZ-QED ALP (left) with tan β = 1 and KSVZ ALP (right).

Figure 7 :
Figure 7: Sketch of the different regions for the angle θ, depending on the ALP-pion mass difference and the sign of ϵ.

Table 1 :
"Experimental" cross sections at BESIII and Belle II for e + e − → V → γa, compared to the non-resonant ones, e + e − → γ * → γa.Vanishing ALP couplings to cand b-quarks have been assumed here.