Revisiting the production of ALPs at B-factories

In this paper, the production of Axion-Like Particles (ALPs) at $B$-factories via the process $e^+e^- \to \gamma a$ is revisited. To this purpose, the relevant cross-section is computed via an effective Lagrangian with simultaneous ALP couplings to $b$-quarks and photons. The interplay between resonant and non-resonant contributions is shown to be relevant for experiments operating at $\sqrt{s}=m_{\Upsilon(nS)}$, with $n=1,2,3$, while the non-resonant one dominates at $\Upsilon(4S)$. These effects imply that the experimental searches performed at different quarkonia resonances are sensitive to complementary combinations of ALP couplings. To illustrate these results, constraints from existing BaBar and Belle data on ALPs decaying into invisible final states are derived, and the prospects for the Belle-II experiment are discussed.


Introduction
Light pseudoscalar particles naturally arise in many extensions of the Standard Model (SM), including the ones endowed with an approximate global symmetry spontaneously broken at a given scale, f a . Sharing a common nature with the QCD axion [1][2][3], (pseudo) Nambu-Goldstone bosons are generically referred to as Axion-Like Particles (ALPs). The ALP mass m a can, in general, be much lighter than the symmetry breaking scale f a , as it is paradigmatically exemplified in the KSVZ and DFSZ invisible axion models [4][5][6][7]. Therefore, it may be not inconceivable that the first hint of new physics at (or above) the TeV scale could be the discovery of a light pseudoscalar state.
The ALP parameter space has been intensively explored in several terrestrial facilities, covering a wide energy range [8][9][10][11][12][13][14], as well as by many astrophysical and cosmological probes [15][16][17]. The synergy of these experimental searches allows to access several orders of magnitude in ALP masses and couplings, cf. e.g. ref. [18] and references therein. While astrophysics and cosmology impose severe constraints on ALPs in the sub-KeV mass range, the most efficient probes of weakly-coupled particles in the MeV-GeV range come from experiments acting on the precision frontier [19]. Fixed-target facilities such as NA62 [20] and the proposed SHiP experiment [21] can be very efficient to constrain long-lived particles. Furthermore, the rich ongoing research program in the B-physics experiments at LHCb [22,23] and the B-factories [24,25] offers several possibilities to probe yet unexplored ALP couplings.
The main goal of this paper is to re-examine existing BaBar and Belle flavor-conserving constraints on ALPs, and to identify the most promising experimental searches to be performed at the forthcoming Belle-II experiment. While there have been several studies discussing signatures of ALPs at B-factories [26][27][28][29], many clarifications are still needed.

JHEP06(2019)091
Firstly, the resonant contributions to the ALP production, via the e + e − → Υ(nS) → aγ process, have been overlooked before. As will be shown, these effects can induce numerically significant corrections to experimental searches performed at √ s = m Υ(nS) , with n = 1, 2, 3. Another improvement provided here concerns the theoretical expression for the Υ → γa branching fraction. Previous studies estimate this quantity by considering either the ALP coupling to b-quarks [30], or to gauge bosons [27]. In this paper, it will be shown that the simultaneous presence of both interactions, as expected in the most general framework, gives rise to new interesting phenomenological features.
In order to assess the limits on ALP couplings, one should specify not only the ALP production mechanism, but also its decay products. In this paper, it will be assumed that the ALP does not decay into visible particles. Such a scenario can be easily achieved by assuming a sufficiently large ALP coupling to a stable dark sector, as motivated by several dark matter models. The conclusions related to ALP production are, however, general and they can also be applied to the reinterpretation of experimental searches with visible decays in the detector, as will be discussed in the following.
The remainder of this paper is organized as follows. In section 2, the effective Lagrangian describing the interactions between the ALP and SM particles up to dimension five is introduced. In section 3, the relevant non-resonant and resonant contributions to the e + e − → γa process are computed. Section 4 is devoted to the classification of experimental searches that can be performed at B-factories. The phenomenological implications of these results are illustrated in section 5, by reinterpreting present constraints on the benchmark scenario of an ALP decaying into invisible particles, and by discussing the corresponding prospects for the Belle-II experiment. Conclusions and final remarks are presented in section 6.

ALP effective Lagrangian
The dimension-five effective Lagrangian describing ALP interactions, above the electroweak symmetry breaking scale, can be generically written as [11] where V µν = 1 2 ε µναβ V αβ , c af f and c aV V denote the ALP couplings to fermions and to the SM gauge bosons, V ∈ {g, B, W }, respectively. The ALP mass m a and the scale f a are assumed to be independent parameters, in contrast to the QCD-axion paradigm, which is characterized by the relation m a f a ≈ m π f π [31]. Moreover, if the ultraviolet completion of eq. (2.1) is not specified, the ALP couplings to gauge bosons and fermions are described by independent parameters, which can be of the same order of magnitude and which should, therefore, be simultaneously considered in phenomenological analyses.
At the energy-scales relevant at B-factories, the ALP interactions with the Z boson can be safely neglected, due to the Fermi constant suppression. Furthermore, the ALP JHEP06(2019)091 couplings to the top-quark and W ± boson are relevant only to the study of flavor-changing neutral currents observables, which are complementary to the probes discussed heresee e.g. refs. [32][33][34][35] for a recent discussion. The only relevant couplings in eq. (2.1) at low-energies are where c aγγ = c aBB cos 2 θ W + c aW W sin 2 θ W . The couplings relevant to ALP production are {c aγγ , c abb }, while the other couplings only contribute to the ALP branching fractions. Light pseudoscalar particles can also act as portals to a light dark sector [36,37]. In this case, to describe these additional interactions, new couplings are customarily introduced. By assuming, for instance, an extra light and neutral dark fermion state χ, the following term should be considered in the effective Lagrangian: where c aχχ denotes a generic coupling, which can induce a sizable ALP decay into invisible final states, as will be considered in the following.
In the remained of this paper, c aii ≡ c eff aii (µ = m b ) will be assumed, and the ALP mass will be taken in the range m a ∈ (0.1-10) GeV, for which B-factories provide some of the most stringent bounds on its couplings.

B-factories probes of invisible ALPs
In this section, the potential of B-factories to probe ALP couplings in the e + e − → γa channel will be discussed. Two main scenarios are typically considered in the literature, depending on the relative strength of the ALP coupling to SM and dark sector particles: either |c aχχ | |c aSM |, or |c aχχ | |c aSM |. In the first case, for m a values in the GeV range, the ALP would typically decay in the detector, leaving the signatures γa(→ jj) [38,39], γa(→ γγ) [40][41][42] and γa(→ ) [43][44][45], with = {e, µ, τ }. This scenario is dubbed the visible ALP. If, however, the coupling to the dark sector c aχχ is large, in comparison to the SM couplings, then the ALP will decay predominantly into an invisible channel, providing the mono-γ plus missing energy signature. This scenario will be referred to as the invisible ALP, 1 which also covers the possibility of a sufficiently long-lived ALP that does not decay in the detector.
In this paper, the invisible ALP scenario will be considered for the sake of illustration. The main goal will be (i) to revisit the theoretical expressions available in the literature, including ALP coupling to bottom quarks, as well as previously unaccounted experimental uncertainties, and (ii) to propose an optimal strategy for future ALP analyses. Even though the main focus will be the minimalistic invisible ALP scenario, most of the observations that will be made in this paper can be translated mutatis mutandis to the visible case. Non-resonant ALP production. The most straightforward way of producing ALPs in e + e − facilities is via the non-resonant process e + e − → γa, as illustrated in figure 1. If the ALP does not decay inside the detector, as assumed throughout the paper, this process would result in an energetic γ plus missing energy. The differential cross-section for this process, keeping explicit the ALP mass dependence, can be expressed as [46] dσ(s) d cos θ γ NR = α em 128 where s = E 2 cm , and θ γ is the angle of photon emission with respect to the collision axis, in the center-of-mass frame. In this expression, the contributions coming from the exchange of an off-shell Z boson, which are also induced by c aW W in eq. (2.1), have been neglected, since they are suppressed, at low-energies, by s/m 2 Z 1. The integrated cross-section then gives: While the non-resonant contribution to ALP production given above is unavoidable in any experiment relying on e + e − collisions [26], the situation at B-factories is more intricate since these experiments operate at specific Υ(nS) resonances. Therefore, it is crucial to account for the resonantly enhanced contributions, which can be numerically significant, as will be shown in the following.
Resonant ALP production. Vector quarkonia can produce significant resonant contributions to the mono-γ channel, e + e − → Υ → γa, since they are very narrow particles coupled to the electromagnetic current. Assuming a fixed center-of-mass energy √ s ≈ m Υ , as is the case at B-factories, and using the Breit-Wigner approximation, one finds for the resonant cross-section where m Υ and Γ Υ are the mass and width of a specific Υ resonance, and σ peak is the peak cross-section defined as with B(Υ → ee) being the leptonic branching fraction, experimentally determined for the different Υ(nS) resonances [47]. The effective couplings defined in eq. (2.2) appear, instead, in the B(Υ → γa) branching fraction, as illustrated in figure 2, which will be computed in full generality in sections 3.1 and 3.2.
Non-resonant vs. resonant ALP production. Naively, one would expect that the resonant cross-section (eq. (3.3)) clearly dominates over the non-resonant one (eq. (3.2)) for the very narrow Υ(1S), Υ(2S) and Υ(3S) resonances, since Γ Υ /m Υ 1. Nevertheless, this turns out to not be the case at B-factories, since these experiments are intrinsically limited by the energy spread of the e + e − beam, which is of order σ W ≈ 5 MeV at current facilities. 2 This value is considerably larger than the width of these resonances, which therefore cannot be fully resolved at B-factories. The only exception is the Υ(4S) resonance, for which Γ Υ(4S) = 20.5 MeV [50]. Therefore, one should expect a sizable reduction of the estimation in eq. (3.3) for the lightest quarkonia resonances, due to this intrinsic experimental uncertainty.
To account for the beam-energy uncertainties in eq. (3.3), the procedure presented in ref. [51] has been adopted by performing a convolution of σ R (s) with a Gaussian distribution, with spread σ W , (3.5) At the very narrow Υ(nS) resonances, with n = 1, 2, 3, one finds Γ Υ σ W , in such a way that the previous expression can be simplified by writing [51] where the parameter ρ, defined as accounts for the cross-section suppression at the peak due to the finite beam-energy spread. These effects will be quantified in the following in two scenarios: (i) ALP with predominant couplings to photons, |c aγγ | |c abb |, and (ii) the general case with both c aγγ and c abb nonzero.

3.1
The photo-philic scenario: |c aγγ | |c abb | The scenario most commonly considered in the literature is the one with predominant ALP couplings to photons [27]. In this case, by neglecting c abb , the first diagram in figure 2 leads to which agrees with ref. [26] in the massless ALP limit. Note that this expression does not require assumptions on hadronic uncertainties, since the hadronic matrix element appearing in this computation, namely 0|bγ µ b|Υ(p) , also enters the process Υ → ee which has been accurately measured experimentally [47]. Alternatively, B(Υ → ee) can be expressed in terms of the Υ decay constant, defined as which encapsulates the QCD dynamics of this process and which has been independently computed, for the lighter Υ resonances, by means of numerical simulations of QCD on the lattice [52][53][54]. The previous definition allows to recast eq. (3.8) in the more convenient form, (3.10) The Lattice QCD (LQCD) determinations of f Υ(nS) are summarized in table 1 along with the values extracted from the experimentally determined B(Υ(nS) → e + e − ) [47], showing a reasonable agreement. In table 2, eq. (3.8) is combined with eqs. (3.2) and (3.6) to estimate the resonant and non-resonant cross-section, for each Υ resonance, along with the peak cross-section σ peak and the suppression parameter ρ. This computation has been performed with the Belle-II (KEKB) energy-spread for illustration, which is similar to the ones from BaBar (PEP) and Belle (KEK). From this table, one learns that even though the peak cross-section is large for the Υ(nS) resonances (n = 1, 2, 3), the beam-energy uncertainties entail a considerable suppression of the visible cross-section. These effects are milder for the Υ(4S) resonance, but in turn the cross-section at the peak is much smaller in this case. The final results are summarized in the last column of table 2, which shows that the effective resonant cross-section is smaller than the non-resonant one, but it still contributes with numerically JHEP06(2019)091   1, 2, 3), the resonant contribution amounts to corrections between 20% and 50% to the non-resonant one, which should be included when reinterpreting experimental searches. 3 On the other hand, for the Υ(4S) resonance the resonant contribution turns out to be negligible, due to its larger width, as expected.
3.2 The general case: c aγγ = 0 and c abb = 0 The previous discussion implies that the resonant contributions are not only important to correctly assess limits on the ALP coupling to photons, c aγγ , but they also open the window to probe the ALP coupling to b-quarks, c abb , cf. figure 2. The simultaneous presence of these contributions gives rise to a rich phenomenology which will be discussed in the following. Firstly, the hadronic matrix element needed to estimate the c abb contribution to B(Υ → γa) is far more intricate than the one given in eq. (3.9), since this is a QCD-structure dependent emission, as depicted in the last two diagrams in figure 2. This contribution was first computed by Wilczek for a SM-like Higgs by using a non-relativistic approximation [2,55], see also refs. [56][57][58][59]. 4 By using a similar approach, the total B(Υ → γa) branching fraction reads (3.11) This expression includes, for the first time, the most general c aγγ and c abb contributions, as well as their interference. Note, however, that the computation of the c abb contributions are done within a first approximation that considerably simplifies the QCD structuredependent emission of this decay. If a new physics signal is indeed observed in such observable, a more accurate theoretical calculation would be needed to fully assess the (non-perturbative) effects associated to the last two diagrams in figure 2.
3 Interference effects between the non-resonant and resonant caγγ terms turn out to be negligible due to the small width of the Υ(nS) resonances. 4 Compatible results have also been obtained in ref. [60] for small pseudoscalar masses by using a QCD sum-rules approach.

JHEP06(2019)091
As shown in eq. (3.11), the c aγγ and c abb couplings can induce comparable contributions to the non-resonant cross-section in eq. (3.3). Moreover, depending on the relative sign of these two couplings, these couplings can interfere destructively or constructively, as will be illustrated with a concrete example in section 5. Finally, note that eq. (3.11) shows a different dependence on m a and {c aγγ , c abb } than the non-resonant cross-section in eq. (3.2). A comparison between σ R vis and σ NR vis ≈ σ NR is postponed to section 5 where a concrete scenario will be considered.

Summary of experimental searches
From the previous discussion, one learns that the non-resonant cross-section, via the coupling c aγγ , is the largest one, but it can be of the same order of the resonant one, cf. table 2. Moreover, the latter searches have the advantage of being sensitive to both c aγγ and c abb couplings. Based on these observations, ALP searches at B-factories can be classified in the following three categories: i) Resonant searches. Excited quarkonia states Υ(nS) (with n > 1) can decay into lighter Υ(nS) resonances via pion emission, as for example Υ(2S) → Υ(1S) π + π − and Υ(3S) → Υ(1S) π + π − . By exploiting the kinematics of these processes one can reconstruct the Υ(1S) meson and then study its decay into a specific final state, which can, for instance, be the invisible Υ decay [61], or the Υ decay into photon and a light (pseudo)scalar particle [40,41]. These searches are dubbed resonant, since they allow to directly probe B(Υ → γa) in a model-independent way, regardless of the non-resonant contribution from figure 1. In other words, reported limits on B(Υ(1S) → γa) can be used to constrain both c aγγ and c abb via eq. (3.11). Searches along these lines have been performed, for instance, by BaBar [40] and, more recently, by Belle [41], under the assumption that the ALP does not decay into visible particles inside the detector.
ii) Mixed (non-)resonant searches. Alternatively, experimental searches could be performed at Υ(nS) (with n = 1, 2, 3) without identifying the Υ decay from a secondary vertex. Example of such experimental searches are the ones performed at √ s = m Υ(3S) [42], where limits on B(Υ(3S) → γa) × B(a → inv) are extracted from the total e + e − → γa(→ inv) cross-section. From the above discussion, however, it is clear that this method is probing both resonant (eq. (3.6)) and non-resonant (eq. (3.2)) cross-sections and therefore model-independent limits on B(Υ(3S) → γa) could not be extracted from these experimental results. The only scenarios for which such limits can be derived are the ones with |c aγγ | |c abb |, as predicted in models with an extended Higgs sector [62][63][64], since the non-resonant cross-section vanishes in this case.
In the most general ALP scenario, instead, the limits on {c aγγ , c abb } can be obtained from ref. [42] via a rescaling factor,

JHEP06(2019)091
which accounts for the non-resonant contributions (eq. (3.5)) that have been overlooked experimentally in the total cross-section. For instance, in the case where c abb = 0, one obtains constraints on c aγγ which are a factor of ≈ 3 more stringent than the estimation that overlooks the latter effects, cf. table 2. Note, also, that similar effects have also been overlooked in reinterpretations of other experimental limits, as for example the ones on B(Υ(3S) → γa) × B(a → hadrons) [38] to constrain the product of ALP couplings to photons and gluons [28].
The reinterpretation described above, for the results from ref. [42] and similar searches, has a possible caveat related to the treatment of the background. In these experimental analyses, the background is determined by considering an independent data sample collected outside the resonance region, typically ≈ 30 MeV below m Υ(3S) . While this strategy allows for a robust determination of the SM background in scenarios with c aγγ = 0, this is not an efficient method if c aγγ is non-negligible. In the latter case, the background sample also receives contributions from the non-resonant diagram in figure 1, which turns out to be the dominant effect. For that reason, it is important that future experimental searches determine the background without relying on off-resonance samples, as performed, for instance, in dark photon searches [65]. Furthermore, it would be helpful to also report the limits on the e + e − → γa crosssection instead of the Υ(nS) branching fraction, as these results contain the full information on both resonant and non-resonant contributions.
iii) Non-resonant searches. The resonant cross-section is negligible at the Υ(4S) resonance, as can be seen from table 2, since its mass lies just above the BB production threshold. Therefore, experimental searches at the Υ(4S) resonance can only probe the c aγγ coupling via the non-resonant ALP production illustrated in figure 1. To our knowledge, no such ALP search has been performed yet at B-factories. For the future, this type of search could exploit the large luminosity collected at Υ(4S) Belle-II, providing the most stringent limits on c aγγ for a GeV ALP mass. See ref. [27] for a recent discussion on Belle-II prospects.
In summary, ALP production receives both resonant and non-resonant contributions at B-factories. The interplay between these production mechanisms allows to classify three complementary experimental strategies: (i) resonant searches of Υ → γa, from which one could infer bounds on c abb and c aγγ , (ii) mixed (non-)resonant searches which are sensitive to a different combination of c abb and c aγγ , and (iii) non-resonant searches which depend solely on c aγγ . Before deriving constraints on the ALP couplings from existing BaBar and Belle data, it is important to stress once again that the conclusions outlined above are general and that they apply, for instance, to searches for ALP decaying into visible particles, such as hadrons [38,39], µµ [43,44] and τ τ [45].

Constraining the ALP parameter space
In this section, constraints on the ALP parameter space are derived from existing BaBar and Belle data, and prospects for the Belle-II experiment are discussed.   [40,42] and Belle [41] are considered. For the constraint on c aγγ obtained from data collected at Υ(3S) resonance [42], the reinterpretation of BaBar limits that neglects the non-resonant ALP production (blue dashed-dotted line) is also considered, along with the rescaled limit that accounts for both resonant and non-resonant ALP production (solid blue line), cf. eq. (4.1). The latter results provide the most stringent limits on the ALP photon coupling from searches at B-factories.
the invisible ALP scenario will be considered, by assuming that B(a → inv) = 1, or equivalently, that the ALP does not decay inside the detector. As anticipated above, the results derived below can be easily recast to scenarios in which the invisible ALP branching fraction is smaller than one.
Firstly, separate constraints on c aγγ and c abb are derived by assuming that the other Wilson coefficient vanishes. These couplings are subject to the limits on B(Υ(1S) → γa) × B(a → inv) reported by BaBar [40] and Belle [41], in which the quarkonia state is produced via the Υ(2S) → Υ(1S)π + π − decay, cf. discussion in section 4. Limits on B(Υ(3S) → γa) × B(a → inv) reported by BaBar [42] are also considered by including the non-resonant contribution overlooked in the experimental analysis, cf. eq. (4.1). These constraints are combined in figure 3 to constrain c aγγ /f a and c abb /f a as a function of the ALP mass. While the limits on c abb from the different experimental searches turn out to be similar, the recast described above of Υ(3S) data provides the most stringent limit on c aγγ . For comparison, the limits obtained by neglecting the non-resonant contribution are also depicted in the same plot by the dashed-dotted line, which turn out to be weaker, as expected. It should be stressed that this reinterpretation is not strictly correct due to the background treatment in ref. [42], but it can be seen as the expected sensitivity of such searches if the background is determined without relying on off-resonance samples, as discussed in section 4.
Next, the allowed parameter space in the plane {c aγγ , c abb }/f a when both couplings are simultaneously considered is shown in  such a way that both couplings interfere destructively in eq. (3.11). In this case, it can be seen from figure 4 that the Υ(1S) constraints have an unconstrained direction that cannot be resolved by only relying on resonant ALP searches. 5 The combination of couplings that lead to this cancellation depends on the ALP mass, especially for m a values near the kinematical threshold, as depicted in the right panel of figure 4. BaBar results obtained at the Υ(3S) resonance, which is not reconstructed, depicts a different sensitivity to {c aγγ , c abb }, as shown by the blue regions in the same plot. While a cancellation between c aγγ and c abb is possible for resonant cross-section, this cannot occur for the non-resonant one (3.2), which depends only on the c aγγ coupling. The combination of these complementary searches allows one to corner the ALP parameter space as depicted in figure 4. Moreover, projections for searches performed at Belle-II, operating at the Υ(4S) resonance, as computed in ref. [27], are displayed in the same plot for an expected luminosity of 20 fb −1 .
Before concluding, comments on studies providing similar constraints on ALP couplings are needed. The authors of ref. [27] have performed a reinterpretation of the BaBar darkphoton search in the e + e − → γ + inv channel [65]. The constraints on c aγγ they obtain, by only considering the non-resonant process from figure 1, are a factor of ≈ 2 better than the limits derived in this paper. Nonetheless, such reinterpretation should be performed with caution for two reasons. Firstly, the kinematical distribution of this process is different for ALPs and dark photons scenarios, as can be inferred from the comparison between eq. (3.2) with the expressions given in ref. [66]. Therefore, to translate the dark photon constraints JHEP06(2019)091 into limits on ALP couplings, one should correct for the different detector efficiencies for the two cases. Another important issue is the fact that the dark photon analysis from ref. [65] combine off-resonance data with data collected at the Υ(2S), Υ(3S) and Υ(4S) resonances. While the photons accompanied by dark photons cannot be produced via Υ(nS) decays, this is not the case for ALPs, as discussed above. Therefore, it is important to account for the resonant ALP production estimated in table 2, which is different for each data set considered by BaBar and which can amount to corrections of O(50%) to the total cross-section.

Conclusions
In this paper, ALP production in association with photons at B-factories is revisited. In particular, the contributions to the e + e − → γa cross section are derived, assuming generic non-vanishing ALP couplings with both photons and b-quarks. The production of ALPs can proceed through the non-resonant channel, e + e − → γa, as well as the resonant one, e + e − → Υ(nS) → γa, which has the unique potential to probe the ALP coupling to b-quarks. After computing the relevant cross-sections and accounting for the effects stemming from the beam-energy uncertainty at B-factories, three distinct and complementary experimental searches have been identified: i) Resonant searches that exploit decays such as Υ(2S) → Υ(1S) π + π − and/or Υ(3S) → Υ(1S) π + π − to directly probe the Υ(1S) decays [40,41], which turn out to be equally sensitive to ALP couplings to photons and bottom quarks, as shown in eq. (3.11); ii) Mixed (non-)resonant searches that use, instead, the primarly produced Υ(nS) resonance, with n = 1, 2, 3, as in the analysis performed in ref. [42]. These searches can probe both resonant and non-resonant ALP production, and hence are more sensitive to the ALP coupling to photons than to the one with b-quarks, cf. section 4; iii) Non-resonant searches, as the ones performed at √ s = m Υ(4S) , that can provide information only on the ALP coupling to photons, cf. table 2. Note, in particular, that neither Babar or Belle have reported such an analysis thus far.
Previous phenomenological analyses overlooked the distinction between these types of experimental searches, which have been clarified in this paper, and the optimal experimental strategies have also been discussed.
To illustrate the phenomenological implications of the effects mentioned above, the scenario with an ALP decaying into invisible final states has been considered. Constraints on the parameter space {m a ; c aγγ , c abb } have been derived from existing BaBar and Belle data, and projections for Belle-II have been discussed. In particular, constraints from resonant searches have a flat direction due to possible cancellations between c aγγ and c abb contributions in B(Υ(1S) → γa). These flat directions, however, can be removed by existing mixed (non-)resonant searches performed at √ s = m Υ(3S) , due to the interplay between resonant and non-resonant contributions described above. In the future, the Belle-II experiment has the great opportunity to perform a first search at √ s = m Υ(4S) , probing solely the JHEP06(2019)091 c aγγ coupling, and providing a complementary piece of information to the aforementioned constraints. Finally, the invisible ALP scenario has been considered in this paper for sake of illustration, but the conclusions derived above also apply to scenarios where the ALP decays into visible particles. The phenomenological study of more general scenarios, including visible ALP decays, as well as experimental signatures with displaced vertices, will be the object of a future work.