Revised constraints and Belle II sensitivity for visible and invisible axion-like particles

Light pseudoscalars interacting pre-dominantly with Standard Model gauge bosons (so-called axion-like particles or ALPs) occur frequently in extensions of the Standard Model. In this work we review and update existing constraints on ALPs in the keV to GeV mass region from colliders, beam dump experiments and astrophysics. We furthermore provide a detailed calculation of the expected sensitivity of Belle II, which can search for visibly and invisibly decaying ALPs, as well as long-lived ALPs. The Belle II sensitivity is found to be substantially better than previously estimated, covering wide ranges of relevant parameter space. In particular, Belle II can explore an interesting class of dark matter models, in which ALPs mediate the interactions between the Standard Model and dark matter. In these models, the relic abundance can be set via resonant freeze-out, leading to a highly predictive scenario consistent with all existing constraints but testable with single-photon searches at Belle II in the near future.


Introduction
Axions and axion-like particles (ALPs) are a generic feature of many extensions of the Standard Model (SM), occurring for example in most solutions of the strong CP problem (see refs. [1,2] for recent work in this direction), in string compactifications [3,4] and in models with broken supersymmetry (the so-called R-axion [5]). Being Pseudo-Goldstone bosons, they can naturally be light and very weakly coupled, thus evading many of the strong constraints on new physics imposed for example by the LHC. In fact, the most promising way to search for ALPs may be the intensity frontier [6,7], i.e. searches with relatively low energy but very large integrated luminosity/intensity. The coming years promise significant progress in this area, driven first of all by the upcoming Belle II experiment [8], but also for example by NA62 [9] and by longer-term projects such as the planned SHiP facility [10]. These searches may soon shed light on the existence and nature of ALPs.
ALPs with masses below the MeV scale can have a wide range of implications for cosmology and astrophysics [11], affecting for example Big Bang Nucleosynthesis (BBN) [12], the Cosmic Microwave Background (CMB) and the evolution of stars. ALPs can constitute cold dark matter (DM) [13] and have also been considered as possible explanations for a range of astrophysical anomalies, such as the over-efficient cooling of certain classes -1 -
Over the past few years an increasing amount of interest has been paid to ALPs in the MeV to GeV range [19][20][21][22][23][24][25][26][27][28]. In this mass range, ALPs are largely irrelevant for astrophysics and cosmology, but they can have a number of interesting implications for particle physics. For example, ALPs have been considered as an explanation for the anomalous magnetic moment of the muon [26,29,30] or for exotic resonances in nuclear transitions [31]. Moreover, it has been pointed out that ALPs may play a crucial role in electroweak symmetry breaking and the solution of the hierarchy problem [32] via the so-called relaxion mechanism [33]. Finally, ALPs offer an interesting possibility to connect the SM to a potential DM particle in such a way that thermal freeze-out can be reconciled with constraints from direct detection experiments [20,34,35].
Given the level of detail and sophistication of many recent analyses, it is rather surprising that many constraints on the ALP parameter space have not been updated for a very long time. In fact, a number of constraints shown in recent studies are taken directly from the very early work by Massó and Toldrà [36,37], even though a wealth of newer experimental data has since become available. Moreover, the early calculations were never meant to be more than order-of-magnitude estimates and are therefore potentially misleading when compared to the projected sensitivities of future experiments. The aim of this paper is to revisit constraints on couplings of ALPs to photons and hypercharge gauge bosons for ALP masses in the MeV to GeV range and to compare existing constraints to the sensitivity of upcoming searches for ALPs. For this purpose, we will discuss bounds from electron-positron colliders, from electron beam dumps and from Supernova (SN) 1987A. 1 While many of these constraints have been investigated in detail in the context of different models, such as hidden photons or light scalars (see e.g. [38][39][40][41][42][43]), we present the first comprehensive reinterpretation of these constraints in the context of ALPs.
A central part of our work is to calculate in detail the projected sensitivity to ALPs of Belle II. While the primary purpose of Belle II is to study the properties of B-mesons, the experiment is in fact ideally suited for a wide range of new-physics searches. Several previous studies have investigated the potential of Belle II to search for invisible Dark Photon decays [40,41,44], which will significantly extend the range of earlier searches at BaBar [45][46][47]. Here we present the first realistic study of experimental backgrounds and detection efficiencies for ALP searches at Belle II. We point out that previous studies have significantly underestimated the sensitivity of Belle II for single-photon searches due to overly conservative background estimates based on an extrapolation of BaBar results. Our study demonstrates that Belle II can already explore new parameter space with early data and in the long run will be highly complementary to future searches for ALPs at the LHC [28] or at SHiP [22,48]. We conclude that Belle II therefore possesses a unique opportunity to discover both visibly and invisibly decaying ALPs.
The latter case is of particular interest in the context of models where the ALP is responsible for mediating the interactions between DM and SM particles. We focus on JHEP12(2017)094 the case that DM couples pre-dominantly to photons. Strong constraints on this scenario from the CMB and indirect detection experiments can be avoided if DM annihilations during freeze-out are resonantly enhanced. We find that in this scenario the observed DM relic abundance can be reproduced in a well-defined region of parameter space, which is presently allowed by all experimental constraints but can be largely probed in the next few years with single-photon searches at Belle II. This paper is structured as follows. Section 2 provides a general review of the effective interactions of ALPs and establishes the notation used in the rest of this paper. The focus of section 3 is on the discussion of existing constraints on the ALP parameter space. The special case of an ALP coupled to DM is discussed in section 4. Finally, we present the projected sensitivity of Belle II in section 5, before concluding in section 6.

Effective interactions of ALPs
This work focuses on the interactions of a pseudoscalar ALP a with SM gauge bosons. Specifically, we consider the Lagrangian where B µν and W i,µν denote the field strength of U(1) Y and SU(2) L , respectively, and we have defined the dual field strength tensors viaB µν = 1 2 µνρσ B ρσ . The parameters m a and f a denote the ALP mass and decay constant, which we assume to be independent parameters.
We emphasize that the Lagrangian that we consider does not include all terms that would be expected to be present in a general effective field theory description of ALPs [25,28]. In particular, we do not consider interactions between ALPs and SM fermions or interactions between ALPs and gluons. This restriction is well-motivated in models where the interactions between ALPs and the SM arise from new heavy fermions that do not carry colour charge. The reason we are interested in ALPs that do not couple to gluons and fermions is that such interactions typically lead to flavour-changing processes (via penguin diagrams or via mixing with the π 0 ), which are tightly constrained by searches for rare decays [20,27]. ALPs coupling dominantly to the gauge bosons of U(1) Y and SU(2) L , on the other hand, are much harder to probe experimentally and require dedicated experimental search strategies. 2 Many of the constraints that we discuss below will remain valid even if the ALP has additional interactions. In fact, additional interactions are expected to increase the ALP production cross section and hence lead to even stronger bounds. Nevertheless, in some cases the presence of additional interactions may actually weaken the constraints. For example, additional decay modes will decrease the ALP lifetime and therefore potentially suppress constraints from experiments searching for long-lived particles. Additional

JHEP12(2017)094
interactions may also lead to the trapping of ALPs in astrophysical objects, weakening the constraints obtained from such systems. The reader should therefore be careful when applying the bounds presented in this work to more complicated models.
After electroweak symmetry breaking, the two terms in eq. (2.1) induce four different interactions between the ALP and SM gauge bosons: where the field strengths and their duals are defined as above. The individual couplings can be calculated in terms of the parameters introduced above. The two couplings of greatest interest for the purpose of this work are where θ W denotes the Weinberg angle.
If c B and c W are independent parameters, so are g aγγ and g aγZ . In particular, for c B ≈ c W one finds g aγZ g aγγ . Nevertheless, as pointed out in ref. [23], there are potentially strong constraints on c W from loop-induced flavour-changing processes like B → Ka. It is therefore particularly interesting to consider the case where c W c B and hence We will refer to the case c B ∼ c W (and hence g aγZ g aγγ ) as photon coupling and to the case c B c W (and hence g aγZ ∼ −g aγγ ) as hypercharge coupling. We emphasize that in both cases ALPs will also couple to pairs of heavy gauge bosons, but we do not discuss the effect of these couplings further.
The aγγ-interaction is of particular importance, as it determines the lifetime τ a of the ALP. The decay width Γ a = τ −1 a is given by It is worth emphasizing that for g aγγ 1 TeV −1 and for m a 1 GeV this decay width is extremely small and hence the ALP decay length can be very large, in particular if the ALPs are produced with significant boost γ a = E a /m a . For a detector of size L D the fraction of ALPs that decay within the detector is given by 3 If the ALPs escape from the detector before decaying, this can prohibit searches for the decay a → γγ. On the other hand, such long decay lengths can facilitate a different kind of search, which focuses on the missing momentum carried away by the invisible ALP. We will discuss existing results and future prospects for both search strategies below. The same interaction can also be responsible for the production of ALPs, for example in e + e − collisions. There are two different production processes of interest: ALP-strahlung (e + e − → γ * → γ + a) and photon fusion (e + e − → e + e − + a), see figure 1. For the former process (and in the limit m a → 0) the differential cross section with respect to the photon angle in the centre-of-mass (CM) frame is given by [30] dσ which has a mild angular dependence and is notably independent of the CM energy √ s for m a √ s. 4 ALP-strahlung therefore typically leads to a photon with sizeable transverse momentum, which is a promising experimental signature.
The cross section for ALP production via photon fusion can be calculated by replacing the colliding particles by their equivalent photon spectra γ(x) and making use of the ALP production cross section from a pair of photons [22]: Unless m a is close to √ s, ALP production via photon fusion typically dominates over ALP-strahlung. However, the ALPs produced in this way are much harder to detect experimentally, as they carry only little energy and therefore decay into relatively soft photons in the laboratory frame. We will return to the experimental feasibility of searches for ALPs produced in photon fusion in section 5.3.
This work focuses on ALPs with mass below 10 GeV, so that the decay a → γZ is forbidden. The aγZ interaction nevertheless plays an important role, as it leads to the decay Z → γ + a [28,49] with partial decay width given by Depending on the ALP lifetime, this process can either lead to the signature Z → γ + inv or to Z → 3γ, both of which can be tightly constrained by experiments.  . Existing constraints on ALPs with photon coupling (left) and hypercharge coupling (right). Proton beam dump constraints are taken from ref. [22], LEP constraints on e + e − → γγ from ref. [21], CDF constraints on Z → γγ from ref. [28], bounds from horizontal branch stars from ref. [11], bounds from visible decays of ALPs produced in SN 1987A from ref [50] and bounds from heavy-ion collisions from ref. [51]. All other constraints have been revisited and updated in the present work.
To conclude this section we note that in principle ALPs may also be produced in Higgs decays, h → Za or h → aa, leading to strong constraints from the non-observation of these decay modes [26,28]. These interactions however only appear when considering effective operators of dimension 6 or higher, and they are not directly linked to the interactions between ALPs and SM gauge bosons. While it is instructive to include these interactions in a general effective field theory approach, they are not generic and may be absent in specific UV completions. We will therefore not consider exotic Higgs decays in this work and instead focus on the phenomenology of the interactions between ALPs and gauge bosons.

Review of bounds on the ALP parameter space
In this section we review existing bounds on the ALP parameter space, updating constraints wherever new data or more precise calculations have become available. Most of the constraints that we will discuss only probe the effective ALP-photon coupling g aγγ . The only exception are constraints from high-energy colliders, which depend on whether the ALP couples to photons or hypercharge. We show a summary of all relevant constraints for both cases in figure 2. All collider and beam dump bounds are provided at 95% confidence level (CL), with the exception of the bounds from BaBar, which are provided at 90% CL. Given that the parameter space under consideration covers many orders of magnitude, the difference between the two choices of CL is imperceptible.

Bounds from electron-positron colliders
Mono-photon searches at LEP. Relevant bounds on the ALP parameter space are obtained from so-called mono-photon searches, i.e. searches for highly-energetic photons in association with missing energy resulting from the process e + e − → γ * → γ + a(→ inv). The bound conventionally shown in the ALP parameter space is often attributed to LEP (see e.g. ref. [19]), but actually goes back to the early analysis from ref. [36] of a monophoton search at the ASP experiment at SLAC [52]. Upon closer inspection, however, it turns out that this search is not sensitive to ALP masses in the sub-GeV region, because it requires E γ < 10 GeV in the fiducial region, whereas for light ALPs the photon energy is given by E γ ≈ √ s/2 = 14.5 GeV. To determine the sensitivity of LEP for ALPs, we follow the re-analysis from ref. [53] of a mono-photon search at DELPHI [54] based on 650 pb −1 at CM energies between 180 GeV and 209 GeV. The best sensitivity for an ALP signal stems from the High Density Projection Chamber, which covers 45 • < θ < 135 • , and from the three highest-energy bins included in the analysis, 0.9 E beam < E γ < 1.05 E beam . We implement detector efficiencies and resolution as detailed in ref. [53] and assume that an ALP will escape unnoticed if it travels more than L D = 260 cm in the radial direction without decaying [55]. 5 We note that mono-photon searches have also been carried out at the LHC (for the most recent analyses see refs. [56,57]), but their sensitivity does not significantly improve on the bound from LEP [19]. Moreover, for these searches the validity of the ALP effective theory becomes a concern [19,25]. We therefore do not show bounds from LHC monophoton searches.
Radiative Upsilon decays. A related class of searches, which can be performed at Bfactories, are searches for radiative decays of Υ(nS) with n = 1, 2, 3. While these searches are typically interpreted in terms of new invisible particles coupling to b-quarks, they also apply to the case where the new particle couples directly to the photon: Υ(nS) → γ * → γ + a. In fact, the corresponding branching ratio is easily calculated in terms of the branching ratio into electrons [36]: where m b denotes the b-quark mass. The bound conventionally shown comes from the Crystal Ball experiment, which gives BR(Υ(1S) → γ + inv) < 4.0 × 10 −5 [58]. However, much stronger bounds can be obtained from more recent measurements, such as the bound BR(Υ(3S) → γ + inv) < 3 × 10 −6 from BaBar [45]. Here we reinterpret this latter constraint under the assumption that the ALP will escape from the detector, if it travels a distance of L D = 275 cm from the interaction point without decaying. 6 5 The event selection includes a veto on energy depositions in any electromagnetic calorimeter. We hence approximate the detector length by the outer radius of the High Density Projection Chamber. 6 The event selection includes a veto of energy depositions in the instrumented flux return (IFR). We assume that a photon has sizeable detection efficiency only if it interacts before it reaches the two outermost IFR absorber layers.

JHEP12(2017)094
Dark Photon searches at BaBar. BaBar has published results from a search for invisibly decaying Dark Photons produced in association with an ordinary photon, e + e − → γA (→ inv) [47]. This search differs from an analogous ALP search in that the photon distribution peaks strongly at small polar angles θ in the CM frame. In the limit m A → 0 one finds [41] dσ where denotes the kinetic mixing parameter. 7 To convert a bound on for Dark Photons into a bound on g aγγ for ALPs we therefore have to correct for the fact that the geometric acceptance will be very different in the two cases. The BaBar analysis considers −0.6 < cos θ < 0.6 for m A > 5.5 GeV and −0.4 < cos θ < 0.6 for m A < 5.5 GeV. By integrating the respective differential cross sections for ALP production and Dark Photon production over these ranges we obtain the fiducial cross section including geometric acceptance. Using these numbers, we can translate bounds on Dark Photons into the ALP parameter space under the assumption that all other selection cuts have the same efficiency for the two models. For very small masses of the invisibly decaying particle, we find that the translation is given by Repeating this calculation for finite ALP masses and taking into account the probability that the ALP decays before leaving the detector (see above) using a detector length of L D = 275 cm [59], we can then reinterpret the full BaBar bound in the context of ALPs.
Radiative Z-boson decays. If g aγZ is non-zero, ALPs can also be produced in the decay Z → γ + a. The resulting experimental signature depends on the ALP decay length. If the decay length is large compared to the size of the detector, the most promising search channel is Z → γ + inv. This process has for example been studied by the L3 collaboration at LEP [60], which quotes an upper limit on the corresponding branching ratio of BR(Z → γ +a) < 1.1×10 −6 . We reinterpret this bound under the assumption that the ALP escapes if it does not decay within L D = 180 cm from the interaction point [61]. 8 Constraints of similar strength are obtained for ALPs decaying close to the interaction point. If the ALP has a high boost but a short lifetime, the two photons produced in its decay will be approximately collinear and may thus appear in the detector as a single photon. The resulting signature can then mimic the forbidden decay Z → 2γ [21,62]. Constraints on this decay mode from LEP have been discussed in ref. [21], but ref. [28] points out that for m a m π 0 even stronger constraints can be derived from the CDF result BR(Z → γγ) < 1.45 × 10 −5 [63]. 7 Note that this expression at face value implies a divergent total cross section. In a more accurate treatment this is regulated by the finite electron mass, which we have neglected. Since realistic detectors will not be sensitive to photon angles close to 0 or π, these details do not matter. 8 The event selection includes a veto of energy depositions in any calorimeter. We hence approximate the detector length by the outer radius of the barrel hadronic calorimeter.

JHEP12(2017)094
For m a 10 GeV, on the other hand, the photons from the ALP decay can easily be distinguished and one obtains the signature Z → 3γ. A recent ATLAS search constrains this decay mode to BR(Z → 3γ) < 2.2 × 10 −6 [64]. Future LHC searches are expected to significantly improve this bound and to extend the sensitivity to lower ALP masses [28].

Bounds from beam dump experiments
ALPs can be produced in electron and proton beam dump experiments via Primakoff production, i.e. the conversion of a photon into an ALP in the vicinity of a nucleus [65]. A number of relevant electron beam dump experiments have been studied in ref. [66], while constraints from proton beam dump experiments have been reviewed in ref. [22]. SLAC E141. Although E141 primarily searched for long-lived particles decaying to e + e − [67], during some of the data taking a photon converter was inserted in front of the detector, giving the experiment sensitivity also for ALPs decaying to photons. The results from this search have been presented in refs. [68,69]. Although results are only provided for a limited range in ALP masses, the given details on the detector geometry are sufficient to extend the search range by rescaling with the appropriate decay probability and assuming that the positron energy spectrum is independent of the ALP mass as long as m a E beam . To perform the rescaling of the decay probability, we need to make an assumption on the typical energy of ALPs that leads to an observable signal, i.e. a positron with E > E beam /2. We find that good agreement with ref. [68] is obtained for E a ≈ 6.5 GeV. Previous attempts to reinterpret data from E141 have assumed that smaller ALP energies are sufficient to produce sufficiently highly-energetic positrons, leading to somewhat more aggressive bounds [6].
SLAC E137. The E137 experiment at SLAC has published results from a dedicated search for ALPs coupling only to photons [70]. The paper does not consider the turnover of the exclusion limit towards large couplings due to the exponential suppression of the number of ALPs that reach the detector, but it is possible to include this additional effect with the information provided in the paper. Specifically, we start from the photon track-length distribution as a function of energy provided in figure 14 of ref. [70] and assume that the cooling water in the beam dump yields the dominant contribution to the ALP production. 9 Although the total decay length between absorber and detector is 204 m, ALPs decaying at the beginning of the decay volume will often not lead to an observable signal, since the photons produced in the decay may miss the detector, which only has a radius r D ≈ 1.5 m. The typical opening angle between the two photons produced in the ALP decay is θ γγ ∼ 2/γ a = 2 m a /E a . If the ALP decay happens at a distance ∆z from the detector, the two photons will thus hit the detector with a typical separation of θ γγ ∆z. To be conservative, we therefore assume that an ALP decay will only lead to an observable JHEP12(2017)094 signal in the detector if ∆z < E a r D /m a , so that the effective length of the decay volume is reduced for ALPs with small boost factor. This procedure allows us to extend the bound from ref [70] to the full parameter space relevant for ALP searches. Our final result roughly agrees with the one shown in ref. [6], although we find slightly stronger constraints towards larger couplings by using a more realistic photon distribution in our calculation.
Proton beam dumps. Constraints on the ALP parameter space from proton beam dumps have recently been studied in ref. [22]. Here we show the bounds obtained in that work from CHARM [71] and NuCal [72,73].

Bounds from heavy-ion collisions
Ref. [24] proposed a novel way to search for ALPs in ultraperipheral heavy-ion collisions. In such collisions ALPs can be produced via the fusion of two coherently emitted photons, i.e. via diagrams analogous to the one shown in the right panel of figure 1. The resulting ALPs have very low boost factors and their decays consequently lead to two very soft photons that are approximately back-to-back. The minimum photon energy required by the trigger places a lower bound on the ALP masses that can be probed with this search strategy. A reinterpretation of the recent ATLAS search for light-by-light scattering [83] leads to relevant constraints for ALPs with m a > 7 GeV [51].

Bounds from astrophysics
Supernova 1987A. Weakly coupled particles such as axions or ALPs with masses up to about 100 MeV can be copiously produced in the hot core of a supernova. Because of their weak couplings these particles stream out of the core and thereby constitute a new energy loss mechanism. In the absence of such new particles the main cooling mechanism is due to neutrino emission. The corresponding neutrino signal has been observed in the case of SN 1987A, placing a bound on possible exotic energy loss mechanisms, which should not exceed the energy loss via neutrino emission.
ALPs that couple exclusively to photons are produced via the Primakoff process. As electrons are highly degenerate in the supernova core, their phase space is Pauli-blocked and their contribution to ALP production is negligible. Protons are only partially degenerate and correspondingly the process γ + p → p + a is the main production mode. The resulting ALP energy spectrum has been calculated with detailed account of the production process in a core-collapse supernova [74]. An analytical fit of these results has been performed in ref. [50], which yields an approximate expression for the ALP production rate: where T is the effective temperature and σ pr (E a ) is the cross section for Primakoff production [11]. For SN 1987A a good fit is obtained for C = 2.54×10 77 MeV −1 and T = 30.6 MeV. The total energy outflow is then given by which should be smaller than ∼ 3 × 10 53 erg. Using this constraint leads to an upper bound on the ALP coupling of around g aγγ < 6 × 10 −9 GeV −1 for small ALP masses. This bound however does not apply for arbitrarily large couplings, because at some point the axions will interact so strongly that they are trapped in the supernova core. The ALP mean free path is given by where n p denotes the proton density and σ bc is the cross section for back-conversion, which is directly related to the production cross section [75]: with β a = 1 − m 2 a /E 2 a . But even if the ALP mean free path is smaller than the size of the supernova core (∼ 10 km), energy transport via ALPs can still be large. The new particle is harmless only if it interacts more strongly than the particles which provide the standard mode of energy transfer, i.e. neutrinos.
To estimate the rate of energy transport, we can calculate the ALP Rosseland mean opacity κ a . For the inverse Primakoff process this can be written as [76] κ while the contribution from decay is given by [77] κ The constraint from SN 1987A only applies if κ a = κ P a + κ D a < κ ν , where κ ν ≈ 8 × 10 −17 cm 2 /g is the neutrino opacity [36]. For small ALP masses this corresponds to g aγγ 7 × 10 −6 GeV −1 .
An additional constraint from SN 1987A can be obtained by considering not only the energy loss due to ALP emission, but also the visible signal resulting from the "ALP burst" if the ALPs subsequently decay into photons [50]. In figure 2 we show the constraint obtained in ref. [50], which extends the bound from SN 1987A to even smaller couplings.
We emphasize that our treatment of the limit from SN1987A remains somewhat simplistic. Our parametrisation of the ALP energy spectrum and outflow are based on a fit to simulations performed with g aγγ = 10 −10 GeV −1 for a light ALP, and it is unclear how accurate this is over for large values g aγγ and for m a near the supernova core temperature. Modern models of supernovae (see ref. [78]) are also considerably more sophisticated than ours, which assumes a fixed nuclear density and temperature (following the treatment in refs. [36,79]). The mean free path of neutrinos depends sensitively on radius, time and energy, making it difficult to encode the effect in terms of a single opacity [80]. We estimate this introduces a theoretical uncertainty of at least a factor of two into our limits. To reduce this uncertainty will likely require a dedicated simulation of heat transfer in the core of a supernova. While some work in this direction has been undertaken for axions [74,81], it would be interesting to perform dedicated simulations for the ALP scenario.

JHEP12(2017)094
Horizontal branch stars. Constraints on the ALP parameter space can be obtained by considering the cooling of horizontal branch (HB) stars [76,79]. These constraints have recently been investigated in ref. [11], and we show the bounds derived there. As with supernovae, it would be interesting to undertake more detailed simulations for ALPs, as has been done for axions in ref. [82].
Big Bang Nucleosynthesis. It was pointed out in ref. [11] that the ALP parameter space is potentially strongly constrained by BBN, if the ALP decay rate is comparable to the Hubble rate during BBN. In this case, ALP decays into photons may alter the baryonto-photon ratio and the number of relativistic degrees of freedom, leading in particular to conflict with the well-measured ratio of the abundances of deuterium and helium. These constraints were updated in ref. [12] and shown to strongly disfavour the triangular region between the bounds from beam dump experiments and the constraints from SN 1987A and HB stars. Nevertheless, these constraints rely on the assumptions of a standard cosmological history and can be significantly weakened for example if additional relativistic degrees of freedom are present during BBN, making them more model-dependent than the other constraints considered here. We will therefore not show the BBN constraints in the following.

ALPs coupled to dark matter
In this section we will extend the model discussed above and consider the case of ALPs coupled to DM particles. The immediate consequence of such a coupling is that, provided the DM particle is sufficiently light, the ALP obtains an invisible decay mode. This will enhance the sensitivity of searches based on missing energy and suppress the sensitivity of searches that rely on the reconstruction of a visible final state.
For concreteness let us assume that the DM particle is a Majorana fermion χ. The generic interaction with ALPs is then of the form where g aχχ has mass dimension −1, just like g aγγ . The invisible decay width is then given by where m χ denotes the DM mass and Γ inv = 0 for m χ > m a /2. If the coupling g aχχ is large compared to g aγγ and m χ is not much lighter than m a , the invisible decay width can easily become large relative to the visible one. 10 Let us assume that this is indeed the case and that BR(a → χχ) ≈ 100%. Relevant experimental constraints then come from mono-photon searches at LEP and BaBar (as well as from LEP bounds on Z → γ + inv if the ALP couples to hypercharge). The

JHEP12(2017)094
constraints from SN 1987A will also become much stronger, as the DM particles produced in ALP decays can contribute to supernova cooling even in parameter regions where the ALP decay length and mean free path are small. 11 Before we compare the sensitivity of the different searches, let us briefly discuss a related question: is it possible within the simple model introduced above to reproduce the observed DM relic abundance while at the same time evading constraints from various DM searches? If the process χχ → γγ has a cross section close to the thermal one, (σv) th ≈ 3 × 10 −26 cm 3 /s, the DM particle can in principle obtain its relic abundance from thermal freeze-out. For m χ m a the annihilation cross section is given by where v denotes the relative velocity of the two DM particles in the CM frame. This cross section is however tiny even for optimistic parameter choices. For example, choosing g aγγ = 10 −3 GeV −1 , g aχχ = 10 −2 GeV −1 , m χ = 1 GeV and m a = 5 GeV, one obtains σv ≈ 2 × 10 −30 cm 3 /s. A second difficulty arises because the annihilation cross section is s-wave (velocity unsuppressed). Consequently, if σv = (σv) th , the annihilation of DM into mono-energetic photons should still be observable in the present Universe. However, searches for γ-ray lines from Fermi-LAT exclude the thermal cross section for s-wave annihilation into photons across the entire energy range [84]. 12 Thus, thermal freeze-out into photons is only viable if the annihilation cross section has a strong velocity dependence.
A simple way to both enhance the annihilation cross section and to introduce a strong velocity dependence is to consider the case that the DM mass is close to resonance m χ ≈ m a /2. In this case, one needs to consider the full expression for the annihilation cross section: where √ s denotes the CM energy. The quantity relevant for the calculation of the relic abundance is then the thermally averaged annihilation cross section [87] σv = where T denotes the temperature of the thermal bath and K i denote the spherical Bessel functions of the second kind. 11 In principle even DM particles could be trapped inside the supernova, if their interaction rate with real photons (χγ → χγ) and with virtual photons (χp → χp + γ) is large enough. The rate of these processes is however suppressed relative to the rate of ALP back-conversion by a factor of g 2 aχχ m 2 χ 1. We find that the DM mean free path is large compared to the size of the supernova as long as mχ gaγγ gaχχ 10 −7 GeV −1 .
12 While the Fermi-LAT collaboration only considers DM masses down to 200 MeV, independent studies have extended these bounds down to 100MeV [85]. We note in passing that constraints from direct detection experiments are strongly suppressed for this scenario [86] and can be neglected.  Figure 3. Comparison of the thermally averaged annihilation cross section σv (orange, solid) and the annihilation cross section in the limit v → 0 (blue, dotted) as a function of the mass ratio r = m χ /m a for m a = 1 GeV and g aγγ = 10 −4 GeV −1 . The red dashed line indicates the simplified expression from eq. (4.6), which is valid close to the resonance. Away from the resonance the cross section depends also on the ALP-DM coupling, which has been set to g aχχ = 10 −3 GeV −1 .
If m χ is close to (but slightly below) m a /2, the integrand in eq. (4.5) is strongly peaked around s ≈ m 2 a , such that the intermediate ALP is on-shell and the annihilation receives a resonant enhancement. We can then make the replacement s → m 2 a everywhere except in the denominator of eq. (4.4), in which case the integration can be performed analytically. If we furthermore substitute Γ inv from eq. (4.2) for Γ a , we obtain the very simple expression σv π g 2 aγγ x K 1 (x/r) 64 where we have introduced the dimensionless temperature ratio x = m χ /T and the dimensionless mass ratio r = m χ /m a . We make two important observations 1. Contrary to naive expectation, σv does not depend on g aχχ . The reason is that increasing g aχχ enhances the production cross section for the ALP in the intermediate state, but at the same time broadens the resonance. These two effects cancel exactly, so that only the coupling to photons enters in the final expression.
2. Even though σv is exponentially sensitive to the mass ratio r (since K 1 (x/r) ∝ exp(−x/r) √ r), it becomes immediately clear that it is not necessary to tune m χ exactly to m a /2. For example, the difference between r = 0.49 and r = 0.499 is only about a factor of 2 (assuming x = 20, which is a typical value of the temperature ratio during freeze-out). This latter point is illustrated in figure 3, which shows in blue (dotted) the annihilation cross section without thermal averaging, in orange (solid) the annihilation cross section with thermal averaging and in red (dashed) the approximate expression given in eq. (4.6).
As we have seen above, the thermally averaged annihilation cross section (and hence the DM relic abundance) depends only on g aγγ and the mass ratio r. If we fix the mass ratio -14 -

JHEP12(2017)094
r we can therefore determine the value of g aγγ that yields the observed relic abundance. 13 It turns out that the resulting values of g aγγ are in an experimentally interesting region: for r = 0.45 (r = 0.49) we find roughly g aγγ ∼ 10 −4 GeV −1 (g aγγ ∼ 10 −5 GeV −1 ). These estimates indicate the regions of parameter space that are interesting for resonant thermal freeze-out. 14 In the following section we will compare these values to the sensitivity that can be achieved by Belle II (see figure 4).

Sensitivity to ALPs of Belle II
The Belle II experiment at the SuperKEKB accelerator is a second generation B-factory and successor of the Belle and BaBar experiments [8]. It is currently under construction and will start data taking in 2018. SuperKEKB is a circular asymmetric e + e − collider with a nominal collision energy of √ s = 10.58 GeV. The design instantaneous luminosity is 8 × 10 35 cm −2 s −1 , which is about 40 times higher than at the predecessor collider KEKB. The Belle II detector is a large-solid-angle magnetic spectrometer. Three sub-detectors are particularly relevant for the ALP searches described in this paper: a 56-layer central drift chamber (CDC) is used for tracking of charged particles and covers a polar angle region of (17-150) • . The electromagnetic calorimeter (ECL) comprising CsI(Tl) crystals with an upgraded waveform sampling readout for beam background suppression covers a polar angle region of (12-155) • and is located inside a superconducting solenoid coil that provides a 1.5 T magnetic field. The ECL has inefficient gaps between the endcaps and the barrel for polar angles between (31.3-32.2) • and (128.7-130.7) • . An iron flux-return is located outside of the magnet coil and is instrumented with resistive plate chambers and plastic scintillators to mainly detect K 0 L mesons, neutrons, and muons (KLM) that covers a polar angle region of (25-155) • .
We study the sensitivity for a dataset corresponding to an integrated luminosity of 20 fb −1 , which is expected to be collected in 2018 without vertex detectors installed. We then scale the expected sensitivity S(g aγγ ) to the planned full integrated luminosity of 50 ab −1 after about 7 years of running, using S(g aγγ ) ∝ 4 √ L. We argue that the expected increase of beam induced background rates at highest luminosity are not relevant for these searches.
In the following we consider ALP decays into DM and into two photons from ALPs produced in ALP-strahlung (a → χχ and a → γγ) and photon-fusion production (a → γγ only). 13 When calculating the relic abundance by solving the Boltzmann equation, a weak additional dependence on the DM mass arises from the fact that the number of relativistic degrees of freedom and hence the expansion rate of the Universe depends slightly on the temperature during freeze-out. The value of gaγγ implied by the observed relic abundance therefore depends slightly on ma even for fixed mass ratio r. We take this dependence into account by solving the Boltzmann equation numerically using micrOmegas v4.2.5 [88]. We note that this calculation assumes that DM couples sufficiently strongly to photons that it remains in kinetic equilibrium during freeze-out [89].
14 If r is extremely close to 0.5, the ALP-photon coupling can even be somewhat smaller than gaγγ = 10 −5 GeV −1 , but in this case indirect detection constraints again become relevant, as resonant enhancement of annihilation processes can also occur in the present Universe (see figure 3).

ALP decays into dark matter
We study decays of ALPs into DM from ALP-strahlung production for ALP masses up to m a = 8.5 GeV.
Signal Monte Carlo events have been generated using MadGraph5 v2.2.2 [90]. We have generated samples using a fixed ALP mass per sample in steps of 0.05 GeV with 10,000 events each, using a branching ratio into DM of BR(a → χχ) = 1.0. The final state consists of a single, highly energetic photon with an energy where √ s = 10.58 GeV is the collision energy. This search is very similar to the search of Dark Photon decays into DM described in ref. [44]. The backgrounds for this search have been found to be due to high cross section QED processes e + e − → e + e − γ(γ) and e + e − → γγ(γ) where all but one photon are undetected. The background composition is a complicated function of detector geometry details that cannot be adequately reproduced without a full Belle II detector simulation. We therefore take the background rates from ref. [44]. It should be noted that the irreducible background from e + e − → ννγ is negligible. We obtain the signal efficiency for ALPs using generator-level Monte Carlo simulations.
We determine the expected 90 % CL upper limit of signal events n s such that the Poisson probability of observing less than n events when expecting n s + n b events is ≤ 0.1, where n is the integer closest to the number of background events n b . Expected upper limits on the coupling g aγγ are summarized as a function of ALP mass m a in figure 4. The much better expected sensitivity compared to BaBar is mainly due to the more homogeneous calorimeter of Belle II. Figure 4 also shows the parameter ranges corresponding to resonant freeze-out. We observe that, if DM annihilation into photons is resonantly enhanced, existing experiments are not yet sensitive to the values of g aγγ implied by the observed DM relic abundance, but Belle II has a unique potential to probe the parameter regions of particular interest. The sensitivity to high mass ALPs is limited by the trigger threshold for a single photon that will be implemented in Belle II. We conservatively assume a trigger energy threshold of 1.8 GeV which limits the search to ALP masses below m a = 8.6 GeV. If the trigger threshold can be lowered to 1.2 GeV, the sensitivity extends to ALP masses up to m a = 9.3 GeV. A higher collision energy close to the Υ(6S) resonance could further extend the sensitivity to about m a = 9.7 GeV for a trigger threshold of 1.2 GeV.
It should be noted that while the dominant physics background for this study comes from e + e − → γγ(γ) events, the largest fraction of the trigger rate for trigger thresholds 1.8 GeV is due to radiative Bhabha events e + e − → e + e − γ(γ) where both tracks are out of the detector acceptance.

ALP decays into two photons
The experimental signature of the decays into two photons is determined by the relation between mass and coupling of the ALP. This relation affects both the decay length of the ALP and the opening angle of the decay photons. It leads to four different experimental signatures (see figure 5): 1. ALPs with a mass of O(GeV) decay promptly, and the opening angle of the decay photons is large enough that both decay photons can be resolved in the Belle II detector (resolved ).

JHEP12(2017)094
3. Even lighter ALPs decay displaced from the interaction point but still inside the Belle II detector. This is a challenging signature that consists of two reconstructed clusters, one of which has a displaced vertex and contains two merged photons. The latter two conditions typically yield a bad quality of the reconstructed photon candidate which is not included in resolved searches with final state photons. There is however enough detector activity in the ECL or KLM that these are vetoed in searches for invisible final states to reduce high rate e + e − → γγ backgrounds.
4. The lifetime of light ALPs with small couplings is large enough that a significant fraction of ALPs decays outside of Belle II. The experimental signature is a single photon final state that looks identical to the ALP decay into DM (invisible).
The detailed Belle II sensitivity for merged and displaced decays depends strongly on the actual performance of the Belle II reconstruction software and beam background levels which are beyond the scope of this paper. It may be possible to search for displaced clusters in the KLM due to its longitudinal resolution of O(10 cm). The resolved region can potentially be extended towards smaller ALP masses if the ECL reconstruction is improved.
We study resolved ALP decays into two photons (e + e − → γ(a → γγ)) over a range of ALP masses between m a = (0.05-9.0) GeV. Signal Monte Carlo events have been generated using MadGraph5 v2.2.2 [90]. We have generated samples using a fixed ALP mass per sample in steps of 0.05 GeV with 10,000 events each using a branching ratio into photons of BR(a → γγ) = 1.0.
We use a simplified geometry description of the Belle II detector to take into account the ECL geometry acceptance. We use the ECL energy resolution [44] that is expected for 10 % of the full luminosity backgrounds. It should be noted that all three photons have rather high energy and the photon energy resolution is expected to not change significantly for even higher beam backgrounds. We use the approximate ECL crystal positions and sizes to estimate the performance in resolving overlapping photon clusters.
We find that the background is dominated by the QED process e + e − → γγγ with three photons in the final state. Background samples are generated using BABAYAGA.NLO [91][92][93]. Additional small backgrounds for small ALP masses may arise from e + e − → γγ with a third photon candidate coming from beam-induced backgrounds, and from e + e − → γγ where one of the photon converts into an electron-positron pair outside of the tracking detectors. The former will be reduced using the very good time resolution O(ns) of the Belle II ECL at high photon energies [94]. To reduce background from pair conversion, one can use the fact that the secondary electron-positron pair splits in the magnetic field and veto events where the polar angle difference between the photons of the lowest invariant mass photon pair is small and the azimuthal angle difference is rather large. We expect that the pair conversion background is only relevant in the ECL backward region where significant material from the CDC readout electronics is placed about 40 cm away from the crystal front.
A further potential background arises from the SM processes e + e − → π 0 γ, e + e − → ηγ and e + e − → η γ [95]. We therefore exclude +50 MeV −75 MeV mass regions around the nominal η and η masses (our analysis is not sensitive to m a ≈ m π ). In the actual analysis a full study of these backgrounds should be included. Finally, we assume that both beam backgrounds and pair conversion backgrounds can be reduced to a negligible level using the selections described above, without significantly affecting the signal selection efficiency.
Our event selection requires three photons with a CM energy E * > 0.25 GeV and a polar angle in the laboratory frame 17 • < θ lab < 150 • . The invariant mass of the two photons from the ALP decay will peak at the ALP mass. We perform the sensitivity study twice, once using all three possible photon pair combinations (high mass selection) and once using only the photon pair combination with the lowest invariant mass (low mass selection). The latter has a smaller signal efficiency especially at higher ALP masses but a lower combinatorial background. For the three photon combination case we select events where the maximum absolute cosine of the three helicity angles is less than 0.9, and for the two photon combination case we keep events where the absolute cosine of the helicity angle is less than 0.6. These selection criteria maximize the ratio of √ S/B, where S is the number of signal events and B is the number of background events, after all other selection criteria have been applied. It should be noted that the helicity selection criteria not only reduce e + e − → γγγ backgrounds, but will also suppress backgrounds from e + e − → γγ combined with a random third photon from beam backgrounds. We require that all three photons are separated by at least 2 ECL crystals in both polar and azimuthal direction. We do not constrain the three photon invariant mass to the collision energy since our MadGraph signal Monte Carlo does not include additional photon radiation whereas the background Monte Carlo does.
We finally select candidates within [−3σ m 2 , +1.5σ m 2 ] around the generated ALP mass, where σ m 2 is the invariant mass resolution of the decay photon pair. For high mass ALPs we select events within [−3σ γ , +1.5σ γ ] around the expected recoil photon energy (see equation (5.1)) instead. The ranges contains about 85 % of the previously selected signal events. The signal efficiency after all selections is flat and about (35-40) % ((50-55) %) for the two photon (three photon) combination (see figure 6). The photon angle separation distance of 2 ECL crystals is a conservative estimate of the Belle II offline reconstruction performance B el le II 3γ (5 0 ab -1 ) Belle II γ + inv (20 fb -1 ) Belle II γ + inv (50 ab -1 ) LHC g aγZ = -2 tan θ W g aγγ Figure 7. Projected Belle II sensitivity (90 % CL) compared to existing constraints on ALPs with photon coupling (left) and hypercharge coupling (right), as well as the projected sensitivities from SHiP [22] and the LHC [28]. and can likely be improved using advanced reconstruction techniques based on Machine Learning methods, and by using shower shape techniques similar to those applied in high energy π 0 reconstruction. We show the efficiency for single ECL crystal difference for comparison as well.
Events from e + e − → γ(a → γγ) are typically triggered by three energy depositions of at least 0.1 GeV in the ECL. Unlike in the Belle II offline reconstruction, the photon reconstruction at trigger level is much simpler and has a worse angular separation power. We expect that a separation of less than 4 ECL crystals will result in merged photon clusters and make this trigger inefficient for ALP masses below about 0.5 GeV. An ideal trigger will require at least two highly energetic ECL clusters and must not satisfy e + e − → e + e − (Bhabha) vetoes. However, any e + e − → γγ veto decision must be delayed to the high level trigger where offline reconstruction is available in order to maintain a high trigger efficiency for low mass ALPs.
We obtain the expected 90 % CL sensitivity as described above. The sensitivity for long-lived ALPs decaying into two photons is determined from the sensitivity of ALP decays into DM, taking into account the reduced efficiency given by eq. (2.6) using a detector length 16 of L D = 300 cm [96]. The projected sensitivities to the coupling g aγγ are summarized as a function of ALP mass m a in figure 7.
We make a number of important observations from figure 7. First of all, we note that for very light ALPs (i.e. m a ∼ 1 MeV) Belle II single-photon searches can push significantly beyond current constraints from beam dump experiments and can potentially explore the triangular region around g aγγ ∼ 10 −5 GeV −1 , which is currently only constrained by model- Belle II searches for three resolved photons can significantly improve over existing bounds from LEP even with early data (20 fb −1 ). With larger data sets, Belle II will be able to probe couplings of the order of g aγγ ∼ 10 −4 GeV −1 over a wide range of ALP masses. Improvements in the Belle II reconstruction software could push the sensitivity for three resolved photons to slightly lower masses and a dedicated search for displaced photons could extend the long-lived search towards higher masses.
Comparing the two panels in figure 7 we note furthermore that there is a remarkable complementarity between Belle II, SHiP and LHC. SHiP will have greatest sensitivity in the parameter region where the ALP decay length is O(1-100) m, which is difficult to explore with Belle II and the LHC. The LHC, on the other hand, is sensitive mostly to the coupling g aγZ , while Belle II and SHiP directly probe the ALP-photon coupling g aγγ . The combination of these experiments will therefore allow to make significant progress in the exploration of the ALP parameter space. Moreover, we can hope to see an ALP signal in more than one experiment, which would potentially enable us to reconstruct its properties and coupling structure.

Photon fusion
So far we have focused on the case that the ALP is produced in association with a highlyenergetic photon, which facilitates an efficient reconstruction of these events. For ALPs produced in photon fusion the situation becomes more complicated, as the transverse momenta of electron and positron after the collision (and hence their polar angle) are too small to be detectable.
Searches for ALPs produced in photon fusion are interesting for two reasons: first, as shown in the left panel of figure 8 the total ALP production cross section from photon fusion significantly exceeds the one from ALP-strahlung (in particular for small ALP masses), so that photon fusion is responsible for the vast majority of ALPs produced at Belle II [23,30]. And second, the production cross section from photon fusion peaks for small ALP momenta, i.e. ALPs will be produced dominantly at rest (see right panel of figure 8). This means that, in contrast to ALP-strahlung, the opening angle between the two photons produced in the ALP decay will typically be large even for low-mass ALPs.

JHEP12(2017)094
The signature in the Belle II detector will consist of two photons with an invariant mass equal to the ALP mass and missing energy along the beam-pipe. The azimuthal angles of the two photons are back-to-back in the CM frame. The Belle II acceptance for ALPs produced in photon fusion is high: for m a = 0.2 GeV (m a = 2.0 GeV) 66 % (89 %) of all ALPs have both decay photons in the ECL acceptance. However, for low mass ALPs the photon energy is small and often below a typical trigger threshold of 100 MeV per ECL cluster. Studies have shown a very large beam-induced background of low energy ECL clusters [44], making the detection of ALPs produced in photon fusion very difficult. While it may be possible to reject a part of these events in the offline reconstruction, the trigger rates are probably too large to sustain.
A possible opportunity might be to consider the case where either the electron or the positron receives sufficiently large transverse momentum to be tagged [23]. However, the cross section is strongly reduced in this case and backgrounds remain rather large compared to the 3γ final state. A detailed study of the potential sensitivity of Belle II in this search channel is beyond the scope of this work.

Conclusions
In this work we have discussed the experimental situation concerning axion-like particles (ALPs) that interact with Standard Model particles dominantly via couplings to photons and electroweak gauge bosons. Reviewing existing constraints we have argued that the bounds from e + -e − colliders conventionally shown for this scenario are outdated. We have updated these constraints by reinterpreting mono-photon searches at LEP and Dark Photon searches at BaBar. We have furthermore investigated the bounds on ALPs from electron beam dumps and from SN 1987A and provided refined estimates of these constraints. We have also discussed the case that ALPs can be produced in the decay Z → aγ, which is relevant if ALPs couple to hypercharge. A summary of existing constraints is shown in figure 2.
A scenario of particular interest is ALPs coupled to a light DM particle, which induces invisible ALP decays and hence removes a number of experimental constraints. In this case DM can pair-annihilate into photons via ALP exchange. We pointed out that if these annihilations are resonantly enhanced in the Early Universe, the DM particle can be a thermal relic with the required abundance and satisfy all experimental and observational constraints. The DM relic abundance then depends on only two parameters (the ALPphoton coupling and the ratio of DM mass to ALP mass), leading to a highly predictive scenario. Depending on how close the mass ratio is to resonance (m χ ∼ m a /2), the observed DM relic abundance implies ALP-photon couplings in the range g aγγ ∼ 10 −4 -10 −5 GeV −1 . These values are consistent with all existing experimental constraints but lie precisely in the parameter region that can be probed with single-photon searches at Belle II.
Our central observation is that existing bounds on ALPs can be significantly improved with single-photon searches and searches for three resolved photons at Belle II. The former type of search probes invisibly decaying ALPs as well as ALPs that escape from the detector before decaying, while the latter search strategy is highly sensitive to ALPs in the GeV region, which decay promptly into a pair of photons. We have shown that the sensitivity of single-photon searches is substantially better than has been estimated previously due to a significant reduction in background compared to BaBar. In combination these search strategies can cover wide ranges of ALP parameter space and explore various models that are of interest for both cosmology and particle physics (see figures 4 and 7). We have found searches with Belle II to be highly complementary to searches with SHiP (which probes different ALP decay lengths) and at the LHC (which probes different coupling structures).
We have pointed out that there are two regimes that require additional studies with full Belle II simulations in order to understand the sensitivity: ALP decays from displaced vertices and prompt decays of highly-boosted ALPs leading to a merging of the two photons. Finding ways to extend the Belle II sensitivity in these regions will be of great interest for future work. Searching for ALPs produced via photon fusion is one potential avenue to make progress, but backgrounds for this search channel may be prohibitive. Future progress in the search for ALPs therefore requires both the implementation of the search strategies described in this work as well as the exploration of innovative approaches that can reach the remaining corners of parameter space. Combining these with a better understanding of astrophysical constraints of ALPs may finally enable us to understand the role that ALPs play in the Universe.