Model-independent bounds on light pseudoscalars from rare B-meson decays

New light pseudoscalars, such as axion-like particles, appear in many well-motivated extensions of the Standard Model and provide an exciting target for present and future experiments. We study the experimental sensitivity for such particles by revising the CHARM exclusion contour, updating bounds from LHCb and presenting prospects for NA62 and SHiP. We first consider a simplified model of a light pseudoscalar $A$ and then propose a model-independent approach applicable to any spin-0 boson light enough to be produced in B-meson decays. As illustration, we provide upper bounds on $\text{BR}(B \to K\,A) \times \text{BR}(A \to \mu^+\mu^-)$ as a function of the boson lifetime and mass for models that satisfy minimal flavour violation. Our results demonstrate the important complementarity between different experiments resulting from their different geometries.


Introduction
The absence of evidence for new physics at the TeV scale so far has led to a renewed attention for much lighter new particles, which may have evaded detection due to their very weak interactions with Standard Model (SM) particles. One of the most attractive realisations of this idea are light pseudoscalars, often called axion-like particles (ALPs), for which the smallness of the mass and the interactions can simultaneously be explained if they arise as Pseudo-Goldstone bosons of a spontaneously broken approximate global symmetry. These particles may couple to SM particles via the so-called axion portal [1][2][3][4] and can be searched for in a wide range of experiments both at the high-energy frontier of particle physics [5][6][7][8][9] and at the so-called intensity frontier [2,[10][11][12][13].
Light pseudoscalars have received particular attention in the context of dark matter model building. The reason is that if the interactions between dark matter and SM particles are mediated by a pseudoscalar exchange particle, one expects a strong suppression of event rates in direct detection experiments, consistent with the non-observation of a dark matter signal in these exper-iments [1,[14][15][16]. 1 There has consequently been substantial effort to construct and study models with pseudoscalar mediators in the context of dark matter indirect detection [20][21][22][23][24], collider searches [25][26][27] and self-interacting dark matter [4,28].
As pointed out in Refs. [3,10,29,30], models with light pseudoscalars face strong constraints from rare decays. The reason is that any pseudoscalar coupling to SM quarks will at the one-loop level induce flavour-changing processes such as b → s A or s → d A, which can be searched for with great sensitivity in a number of experiments [32]. In the mass range where these constraints are relevant (m A 5 GeV) it is therefore very difficult to obtain sizeable interactions between dark matter and SM quarks via the exchange of a pseudoscalar mediator. Nevertheless, Ref. [30] also found that certain parameter regions are rather difficult to probe experimentally. For example, pseudoscalars with masses between a few hundred MeV and a few GeV and lifetimes of the order of 10 −9 s are very difficult to constrain, since the resulting decay length is of the order of meters, which is large compared to the scales of colliders but short on the scale of fixed-target experiments. It is therefore an interesting and important question to understand which experiments are capable of probing this parameter region. Among the 1 It was recently pointed out that scattering in direct detection experiments may in fact be dominated by loop processes that are unsuppressed in the non-relativistic limit [17][18][19]. But even when including these effects the expected event rates are still so small that a detection will be very challenging. possible contenders are a number of existing and planned fixedtarget experiments with relatively short absorber length, such as NA62 [33,34], SeaQuest [35,36] or SHiP [37,38], collider experiments with sensitivity to displaced vertices, such as LHCb using the vertex locator VELO [39], or dedicated forward detectors like FASER [40,41].
In the present work, we address this question by improving upon the analysis presented in Ref. [30] in four key regards. First, we perform a more careful calculation of constraints and projected sensitivities from fixed-target experiments, taking into account the detailed energy distribution of B-mesons produced in the target as well as the propagation of the resulting particles and the detector geometry using a toy Monte Carlo (MC). Second, we include new constraints from searches for μ + μ − pairs originating from a displaced vertex in rare B-meson decays at LHCb [42,43]. Third, we include updated calculations of the branching ratios of light pseudoscalars from Ref. [44]. Finally, we develop a new approach for presenting experimental bounds that does not require any assumptions on the underlying theory. The aim of all of these improvements is to study in a model-independent way which regions or parameter space remain unexplored and whether they can be probed with future fixed-target experiments.
This work is structured as follows. In Section 2 we introduce the effective interactions that we want to investigate and present the relevant expressions for the partial decay widths of rare Bmeson decays and pseudoscalar decays. Section 3 then concerns the calculation of updated experimental constraints and projections. Our new approach for presenting experimental bounds will be discussed in Section 4, followed by the conclusions in Section 5.

Light pseudoscalars and rare decays
We consider a light pseudoscalar A with Yukawa-like couplings to SM fermions: where m f is the fermion mass, v 246 GeV is the electroweak vacuum expectation value and q = {u, d, s, c, b, t} and = {e, μ, τ } denote SM quarks and charged leptons, respectively. As pointed out by a number of recent studies [21,[25][26][27], the Lagrangian in eq. (1) does not respect the gauge symmetries of the full unbroken SM gauge group and should therefore be thought of as the low-energy limit of a more complete theory, such as a Two-Higgs Doublet Model (2HDM). These theories in general allow for the couplings to down-type quarks to be enhanced by a factor of tan β, but we assume for simplicity that A couples to all SM quarks proportional to their mass. In addition to the flavour-conserving interactions in eq. (1), the electroweak charged current induces flavour-changing couplings at the one-loop level. For the present work we are particularly interested in the one leading to the transition b → s A, which can be written as: with q R,L = 1 2 (1 ± γ 5 )q. As expected for an effective theory, the relevant loop diagrams are divergent, such that the coefficients can be parametrised in terms of a new-physics scale [30]: where m W is the W -boson mass, α ≡ e 2 /(4π ), θ W is the Weinberg angle and V is the CKM matrix. The corresponding expression for h L sb (up to an overall sign) is obtained by replacing m b by m s .
In extensions of the SM there may be additional sources of flavour violation, which can also contribute to h R,L sb . Provided these sources satisfy the assumption of minimal flavour violation [31], the flavour structure will take the same form as eq. (3), but the first factor and the logarithm will differ. We will consider this more general case in Section 4 and construct model-independent ratios of observables. Defining the partial decay widths for the corresponding flavour-changing Bmeson decays are given by [29,32,45,46]: The relevant form factors are given by [47][48][49] D-meson decays involving A, on the other hand, suffer from CKM as well as mass suppression in eq. (3) and are therefore not taken into consideration [50]. Finally, we need to consider the decay of A into SM particles. If the pseudoscalar has no interactions apart from the ones given in eq. (1), the dominant decay modes are expected to be A → + − for most of the mass range that we will be interested in. 2 The corresponding partial decay width is given by where z = m 2 A /(4 m 2 ). The precise branching ratios for light pseudoscalars are however difficult to calculate due to large uncertainties in the partial width for hadronic decay modes [8]. Here we adopt the recent results from Ref. [44], which estimates the hadronic decay width using Chiral Perturbation Theory for m A < 1.2 GeV and a spectator quark model for 1.2 GeV < m A < 3 GeV. 3 For m A > 3 GeV we extrapolate the results from Ref. [44] by assuming that the hadronic decay width converges to the partonic width for a pair of free strange quarks with mass m s = 95 MeV and that decays involving charm quarks can be neglected for 2 In particular, we assume that the pseudoscalar has no enhanced coupling to photons. For a detailed discussion of constraints on ALPs that couple dominantly to photons, we refer to [34,51,52]. 3 Note that, although these results have been obtained in the context of the NMSSM, they can directly be applied to our case. In particular, since the most important contribution to the hadronic decay width stems from the pseudoscalar coupling to strange quarks, the leptonic branching ratios are to first approximation independent of tan β.
GeV. This approach neglects the potential effects of charmonium resonances, which are beyond the scope of the present work. We will return to the question of how to deal with the uncertainties related to the muonic branching ratio in Section 4.

Experimental constraints
In the following we will focus on pseudoscalars that decay predominantly into pairs of muons and taus, i.e. m A 210 MeV, but are light enough to be produced in B-meson decays, m A 4.7 GeV. The dominant constraints in this mass range come from LHCb, which is most sensitive for relatively short-lived pseudoscalars, and fixed-target experiments, which can probe pseudoscalars with very small couplings and rather long lifetime.

Fixed-target experiments
To calculate constraints on light pseudoscalars from fixed-target experiments, the first challenge is the production of B-meson spectra. To this purpose we employ PYTHIA 8.23 [55], solely allowing the simulation of hard-QCD bottom-quark production processes. The PDF set number 2, corresponding to the leading-order 5L set from the CTEQ collaboration, has been used [56]. This number seems in qualitative agreement with the experimental data for higher beam energies shown in figure 30 of [57]. We have cross-checked our spectra against the primary yield computed for SHiP [53] and include with our work an auxiliary file with the spectra. 4 Fig. 1 shows a comparison of the B-meson spectra used in our simulations with the ones found in Refs. [53,54]. As outlined in Ref. [53], in principle it would be possible to also include the yield of B-mesons from secondary proton interaction, which can give a relevant contribution to the overall yield. This contribution is however not included in our study, as it would require much more computationally expensive simulations.
To translate the yield obtained from PYTHIA into a total yield, one has to weigh the MC events with the ratio of the bb cross section and the total proton-proton cross section at 400 GeV. The total proton-proton cross section is taken to be σ (pp) 40 mb, see e.g. Ref. [53] (figure 1). The dependence on the target material is accounted for by a correction factor of A 1/3 [57].
In this study we consider existing bounds from the past CHARM neutrino experiment and prospects for the running NA62 experiment as well as the proposed SHiP facility. 5 To estimate the experimental sensitivities we rely on a toy MC with geometries as outlined below, do not account for efficiencies (unless stated otherwise) and assume that backgrounds are negligible. Our results should therefore not be interpreted as an accurate prediction of what may ultimately be achievable in these experiments after accounting for selection efficiencies and potential backgrounds. 4 The file "B400GeV.root" contains a ROOT Tree of the momenta of B 0 , B 0 S , and B + and their anti-particles, produced after 400 GeV proton beam interactions onto a fixed proton target. The z component is along the beam line. 5 Another fixed-target experiment of interest is SeaQuest [35,36], which operates with a beam energy of 120 GeV. However, to the best of our knowledge, there are no measurements of the bb cross section at a centre-of-mass energy of 15.1 GeV, making it impossible to validate the bb cross section calculated by PYTHIA. We therefore do not include SeaQuest in our analysis. The toy MC for all fixed-target set-ups under consideration has the following structure: First we randomly sample B-mesons (B 0 or B + and their anti-particles) from the PYTHIA output. Then we randomly sample from the decays B → K + A and B → K + A by summing their respective decay widths and accounting for their relative probability. Finally, we calculate the probability that the pseudoscalar decays inside the decay volume of the experiment: where l A = E A τ A /m A is the pseudoscalar decay length, d denotes the distance of the pseudoscalar production point to the decay volume and l is the fiducial length in which pseudoscalar decays can be detected. For CHARM [58,59], the following input enters our toy MC: The The NA62 experiment [60] aims at a precise measurement of the ultra-rare decay K + → π + νν. Besides this main measurement, NA62 has a rich programme for exotic particles, including long-lived particles that can be produced in the up-stream copper beam collimator, where the primary SPS proton-beam is fully (in dump-mode) or partially (in standard data-taking) dumped. Here we consider proton interactions that yield B-mesons in the copperbeam collimator located 23 m downstream the primary target. 6 We model NA62 in our MC according to the following parameters: The distance between the beam-defining collimator and the start of the fiducial volume is d = 82 m and the vacuum decay region before the spectrometer is l = 75 m long. In addition, we require the following acceptance conditions: Both muons produced in the pseudoscalar decay need to be detected at the first and last spectrometer chamber as well as at the calorimeter (for all of these, the central hole that allows the beam to pass has been accounted for). Also, we require that each μ has an energy greater than 5 GeV (in order to be efficiently tracked by the trackers with the standard pattern recognition). Finally, each μ, when extrapolated to the 'MUV3-scintillator' plane, is within its acceptance. The target material for the production of B-mesons is copper and we show NA62 prospects for 1 × 10 18 POT. Finally, we model the prospects for SHiP [37,38] as follows. We implement a spectrometer with the following geometrical parameters [61]: The fiducial region is taken to be 45 m downstream from the production point. The spectrometer magnet is positioned at 32.5 m behind the start of the fiducial region and gives a kick with 0.75 T m. As a minimal requirement, we ask both μ tracks to be within two spectrometer chambers located at 28 m and 37 m behind the start of the fiducial region, respectively. The spectrometer chambers are ellipsoidal regions with half-axes x = 2.5 m and y = 5 m. In addition, both muons need to be detected at the end of a fiducial region of 50 m in the same ellipsoid and we ask for the total energy in acceptance to be greater than 2 GeV. The target material for the production of B-mesons is taken to be molybdenum and we show SHiP prospects for 1 × 10 20 POT.

LHCb
LHCb has performed searches for displaced vertices in the processes B → K φ and B → K * φ, where φ is a scalar boson that subsequently decays into a di-muon pair [42,43]. Since both the B-meson and the scalar boson φ decay isotropically, these searches can be directly reinterpreted for the case of a pseudoscalar instead of a scalar. The results from LHCb then provide upper bounds For B → K * φ LHCb provides a numerical code with an interpolation of the experimental results [62]. For given m A and g Y we therefore simply need to calculate the model prediction for BR(B → K * A) × BR( A → μ + μ − ) and the lifetime τ A and then apply the code to determine whether the point is excluded by LHCb data. For B → K φ, we make use of a digitised version of figure 4 from Ref. [43] to implement a similar approach.
Since events with a di-muon invariant mass in the J /ψ , ψ(2S) and ψ(3770) resonance regions are not considered in the LHCb 6 The '0 m' position at NA62 is the location of a Beryllium target, which is the source of the secondary kaons used for NA62's main measurement. Since NA62 is a 'kaon factory', the consideration of pseudoscalar production in upstream or downstream kaon decays is an interesting topic that goes beyond the scope of this work. analyses, we do not obtain any constraints for pseudoscalar masses close to these resonances. 7 The strongest constraints in these regions therefore still come from the analysis in Ref. [30] of prompt decays in LHCb [63].

Results
The results of our analysis are summarised in Fig. 2. In the upper panel we compare the new exclusion limits (coloured regions) to the ones obtained in previous analyses of light pseudoscalars (grey hatched regions), which only included prompt decays for LHCb and assumed a fixed B-meson energy for CHARM [30]. 8 The shift of the updated CHARM contour towards larger couplings is a direct consequence of the improved calculation of the B-meson energy spectrum. In fact, we find that the B-meson energy for 7 In the analysis of B → K φ, LHCb has also fully vetoed the K 0 S resonance. Moreover, the ψ(4160) resonance is also excluded in the digitised data that we have used. However, in both of these regions we obtain bounds from the B → K * φ data. 8 Note that the definition of g Y adopted in the present work differs by a factor of √ 2 from the one in Ref. [30]. Also, the CHARM bound shown in Ref. [30] included a contribution from the process K → π A, even though kaons loose a substantial fraction of their energy in the target before decaying and therefore do not produce a collimated beam of pseudoscalars. The grey region shown in Fig. 2, in contrast, includes only pseudoscalars produced via B → K A and B → K * A. pseudoscalars in the signal region peaks significantly above 25 GeV, depending on the value of m A . The increased boost factor means that pseudoscalars with larger couplings (and hence shorter lifetime) can still induce a signal in the detector behind the absorber.
Conversely, by including searches for displaced vertices the exclusion limits from LHCb are extended towards significantly longer pseudoscalar lifetimes, i.e. smaller couplings. We find that the improved treatment of CHARM in combination with the updated constraints from LHCb almost completely closes the gap between the two experiments that previously existed around couplings of order 10 −4 . Indeed, for 2m μ < m A 3 GeV the combined constraints imply g Y 5 × 10 −5 .
We also observe a suppression of experimental sensitivities for m A ≈ 958 MeV. This suppression is the result of mixing between the new pseudoscalar and the η -meson, which substantially increases the hadronic decay width of the pseudoscalar. As a result, the muonic branching ratio is suppressed and the pseudoscalar lifetime is reduced, which in combination leads to a suppression of the expected number of observable decays. 9 A similar effect occurs for m A ≈ m η = 548 MeV, but the width of the resulting dip is too narrow to be resolved.
In the lower panel we compare the existing constraints from CHARM and LHCb with the projected sensitivities of NA62 and SHiP. As expected, SHiP promises a substantial gain in sensitivity compared to existing experiments. 10 On the other hand, NA62 on first sight appears not to provide any new information, as the parameter region that can be probed falls exactly into the region excluded by the improved analyses of CHARM and LHCb. Nevertheless, we emphasise that this conclusion only holds in the context of the specific model that we have considered so far. In the following section, we will show by adopting a model-independent approach that NA62 can in fact probe parameter regions excluded by neither CHARM nor LHCb.

Model-independent bounds
The results presented in the previous section have been obtained under a number of assumptions. First of all, we calculated the total decay width of A by adopting a specific calculation of the hadronic decay width and neglecting final states involving charm quarks. The corresponding uncertainties, in particular concerning the mixing between A and pseudoscalar mesons, are however difficult to quantify. Moreover, it is conceivable that the pseudoscalar couples with different strength to quarks and leptons, for example in lepton-specific 2HDMs [64], which would modify our results. An easy way to address these issues and to allow for a broader reinterpretation of our results would be to quote the largest value of BR( A → μ + μ − ) that can be excluded experimentally for each point in g Y -m A parameter space.
We have, however, also made another more troublesome assumption by fixing the relation between the production rate of A and its subsequent lifetime via eq. (3) with a specific choice of . In fact, these two quantities enter into the calculation of experimental bounds in very different ways: The expected number of 9 If we were to include hadronic final states in our analysis, experimental sensitivity might even increase for m A ≈ m η and very small couplings, but a reliable calculation of the hadronic partial width is very difficult in this regime. Lacking an appropriate MC simulation for hadronic final states, we include only muons and thus provide a conservative sensitivity estimate. Note also that we do not include additional contributions to the pseudoscalar production rate resulting from mixing (see also Section 4). 10 In analogy to the improvement for CHARM, we find that the expected sensitivity of SHiP significantly exceeds the one previously published in Ref. [37], where once again a fixed B-meson energy was assumed. events in a given experiment typically depends linearly on the production rate (i.e. the B-meson branching ratios) but exponentially on the pseudoscalar lifetime. It is thus impossible to use the results presented above to infer the corresponding constraints for a different value of using a simple rescaling. At first sight, this problem is not very severe, given that only enters logarithmically. However, it has been shown that the flavour-changing interactions can in fact vary quite substantially between different high-energy theories that lead to the same effective interaction between pseudoscalars and quarks [46]. For example, an important contribution to flavour-changing processes may arise from pseudoscalar-Higgs interactions or interactions between the pseudoscalar and SU (2) L gauge bosons [65]. Furthermore, in specific UV completions there may be additional non-logarithmic contributions, which become important if is small, in particular if they have the opposite sign [66]. Finally, mixing between the pseudoscalar and QCD resonances may enhance the rate of rare B-meson decays for specific pseudoscalar masses.
It therefore becomes an important problem to present results in such a way that they do not rely on a specific relation between production and decay. Indeed, such model-independent approaches are frequently employed in the context of searches for long-lived particles (see e.g. Refs. [67,68]). The crucial idea is to treat the production mechanism and the decay length as fully independent parameters, rather than attempting to derive them from the same fundamental interactions. In our case, this means that rather than presenting experimental bounds as a function of g Y and m A , we present them as a function of μ A ≡ BR(B → K A) × BR( A → μ + μ − ), τ A and m A . The fact that we now have three rather than two independent parameters is a direct consequence of the need to capture a broader range of possible underlying theories.
For a single experiment, the most convenient way to present experimental bounds is to provide upper bounds on μ A as a function of τ A and m A , as done for example in [42,43]. This approach makes it however difficult to compare the sensitivity of different experiments. We therefore prefer to provide upper bounds on μ A as a function of τ A for fixed m A . In principle, this must be done independently for both μ A = BR(B → K A) × BR( A → μ + μ − ) and μ * A = BR(B → K * A) × BR( A → μ + μ − ), since the two B-meson branching ratios depend in different ways on the flavour-changing coefficients defined in eq. (2). Nevertheless, in the case of minimal flavour violation the flavourchanging coefficients must have the same flavour structure as in eq. (3), such that the ratio , (11) does not contain any model-dependent quantities apart from m A . For all models that satisfy minimal flavour violation, we can thus combine constraints on both μ A and μ * A into a single plot. For the case of fixed-target experiments this is done by including both decay modes in the signal simulation, fixing the ratio between them according to eq. (11). In the case of LHCb, on the other hand, we evaluate the constraints on μ A and μ * A separately and show the stronger one.
We present our results in Fig. 3, where the four different panels correspond to different values of m A . As in the bottom panel of Fig. 2, shaded regions correspond to existing exclusion bounds, while dashed and dotted regions indicate the sensitivity of planned experiments. We observe that the different experiments achieve their maximum sensitivity (i.e. the strongest bound on μ A ) for different values of the lifetime τ A , reflecting the different geometries and different Lorentz-boost factors of the pseudoscalars. Indeed, the sensitivity is typically maximised if the pseudoscalar decay length corresponds roughly to the distance of the sensitive volume of the experiment from the production point. Hence, LHCb probes mostly short-lived pseudoscalars, whereas fixed-target experiments achieve the best sensitivity for longer lifetimes.
In particular, we find that the sensitivity of NA62 and CHARM peak at different lifetimes. This is because the B-meson energies, and hence the pseudoscalar Lorentz boosts, are comparable in both experiments but the distance of the sensitive volume is much shorter in NA62 than in CHARM. We see that this different geometry enables NA62 to probe low-mass pseudoscalars with lifetimes 100 ps ≤ μ A ≤ 1000 ps, a region for which neither LHCb nor CHARM are very sensitive. This finding is in stark contrast with the conclusions of Fig. 2, where we found NA62 unable to provide relevant constraints. In other words, the model-independent analysis reveals the unique potential of NA62 to probe certain regions of parameter space. . Consistent with our observations there we find that with this assumption NA62 does not probe any parameter combinations that are not already excluded by the combination of LHCb and CHARM. This conclusion also does not change when changing the value of assumed in the calculation of B-meson decays by moderate amounts. Increasing simply shifts the line towards larger values of μ A , such that the constraints from LHCb and CHARM become even stronger (and vice versa). Adding an invisible decay channel would shift the line to the bottom-left (as both the lifetime and the leptonic branching ratios are decreased), such that experimental constraints can be evaded.
A key advantage of the model-independent approach is that we did not need to assume at any point that the light boson is a pseudoscalar (i.e. CP-odd). Indeed, the constraints that we show apply equally to light scalars that satisfy minimal flavour violation (see e.g. Refs. [69][70][71]). For illustration, we also indicate in Fig. 3 the model-specific predictions for a light scalar that mixes with the Higgs boson (studied most recently in Ref. [69]). This model has the advantage that the loop-induced flavour-changing coefficients are finite and hence independent of the new-physics scale . The calculation of the scalar lifetime and branching ratios, on the other hand, is much more challenging, as hadronic decay modes cannot be neglected. The combination of these two effects shifts the line for the scalar case to the bottom and to the left. Indeed, we find that NA62 may possess a unique sensitivity to this model for scalar lifetimes in the range 100-1000 ps.

Conclusions
Light pseudoscalars arise naturally in many extensions of the Standard Model and may be an important ingredient of phenomenologically viable dark matter models. Important constraints on these models arise from searches for long-lived particles in fixed-target experiments and from searches for rare B-meson decays at LHCb. In this work we have improved the analysis of existing constraints from these experiments and proposed a new method for comparing the sensitivities of future searches.
For the first part of the paper we have adopted an effective low-energy description of the interactions of pseudoscalars with quarks and charged leptons. Loop-induced flavour-changing processes then induce a sensitivity to new physics at higher scales, parametrised by the unknown scale , which enters in the production rate of light pseudoscalars via rare B-meson decays. Evaluating experimental constraints for fixed values of , we find that searches for displaced vertices at LHCb and an improved calculation of the constraints from CHARM in combination rule out much of the relevant parameter space.
Nevertheless, the unknown relation between the rates of production and decay of light pseudoscalars (together with additional theoretical uncertainties in the signal calculation) makes it desirable to adopt a more model-independent approach. To this end we propose to express experimental results and sensitivities as bounds on products of branching ratios like μ A = BR(B → K A) × BR( A → μ + μ − ) as a function of the pseudoscalar lifetime τ A and mass m A . This form of presentation has the additional advantage that it makes the role of different experimental geometries explicit, because the sensitivity is typically maximised for pseudoscalar decay lengths comparable to the size of the experiment.
Adopting this model-independent approach, we demonstrate that NA62 has potentially world-leading sensitivity for pseudoscalar lifetimes in the range of 100-1000 ps. Further substantial progress can be expected from dedicated facilities to search for long-lived particles like SHiP. Our approach can easily be generalised to include other pseudoscalar decay modes, such as A → γ γ or A → 3π , which yield complementary information. In summary, while there exist already strong constraints on light pseudoscalars, substantial terrain remains unexplored and provides room for further exploration and potential discoveries.