Light in the beam dump - ALP production from decay photons in proton beam-dumps

The exploration of long-lived particles in the MeV-GeV region is a formidable task but it may provide us a unique access to dark sectors. Fixed-target facilities with sufficiently energetic and intense proton beams are an ideal tool for this challenge. In this work we show that the production rate of Axion-Like-Particles (ALPs) coupled pre-dominantly to photons receives a significant contribution from daughter-photons of secondary $\pi^0$ and $\eta$ mesons created in the proton shower. We carefully compare the PYTHIA simulated spectra of such secondaries to experimental literature, compute the ALP flux from the Primakoff conversion of these photons, and finally revisit existing limits on ALPs and update the prospects for a set of existing and future searches. Our results show that taking this production mechanism into account significantly enhances the sensitivity compared to previous studies based on coherent ALP production in primary proton-nucleus interactions.


Introduction
At present, the search for new particles is driven by a number of pertinent observations that cannot be explained within the Standard Model (SM). Perhaps most prominent amongst them is Dark Matter (DM), but baryogenesis as well as the strong CP problem also ask for a solution. At the same time, we are challenged by the (so-far) non-detection of new particles. This could be explained by the new particles being extremely weakly coupled. If this is the case they can escape constraints from searches with colliders. But their small coupling also makes such particles long-lived.
ALPs have recently received a considerable amount of interest in the context of dark matter model-building. They may act as a mediator for the interactions between DM and SM particles and thereby allow reproducing the correct Dark Matter relic abundance via thermal freeze-out. At the same time this helps evading the strong constraints from direct and indirect detection experiments [3][4][5]7]. Other motivations for ALPs with masses above O(1) MeV include their potential connection to explaining the observed value of the magnetic moment of the muon [16]. Also, via the so-called 'relaxion' mechanism, ALPs may play a crucial role in electroweak symmetry breaking [17] and in the solution of the hierarchy problem. A concrete implementation [18] of such a relaxion may yield signatures that are observable with the methods described in this paper. Moreover, it is worth reiterating that recently models have been proposed that allow the QCD axion to live at the MeV-GeV scale [12][13][14]. Beyond these phenomenological considerations, ALPs are also motivated by top-down extensions of the SM such as string theory [19][20][21][22]. The crucial feature is that a weakly broken or anomalous shift symmetry allows their mass to be much smaller than the fundamental scale, making them accessible to experimental tests. The same shift symmetry then also ensures that the interaction is suppressed by the fundamental scale and therefore very weak. In summary, new searches for ALPs are a well-motivated and timely task.
The main aim of the present paper is to study the production of ALPs in proton beam dumps from decay photons of secondary mesons. We therefore focus on pseudoscalar ALPs whose dominant interaction is with photons. We use the Lagrangian where g aγ is the photon-ALP coupling and F µν is the electromagnetic field strength.
To discover ALPs with masses on the order of MeV-GeV, proton-fixed target facilities (or rather proton beam-dump experiments) are well suited 1 . The strength of such an experimental setup is that it can provide sufficient energy to produce MeV-GeV scale particles, while ensuring that all of the protons in the beam ultimately interact. Moreover, decay volumes spanning tens of meters allow ALPs of various lifetimes to be detected. Overall this combination provides high sensitivity but also excellent complementarity to experiments at low-energy colliders, such as Belle-II [7], which can explore the region below a few GeV but relatively strong coupling, as well as to experiments at the LHC, which are sensitive mostly to masses above ∼ 5 GeV [28][29][30][31].
For ALPs coupled predominantly to photons, there are (at least) two important production mechanisms to be considered. In both cases, one of the photons is provided by a nucleus at rest, constituting the target/dump material. The second photon can be provided 1. by the charged proton itself. This is sometimes referred to as 'photon-from-proton' (PFP) mode, where the photon distribution around the proton is often computed in the Weizsaecker-Williams approximation 2 2. through a decay, notably from π 0 → γγ, but also from other neutral mesons.
As already noted, the second process is the main focus of our paper. 1 Indeed there is a sizable number of such beams around the world, cf. e.g. [23][24][25][26] whose suitability for ALP searches has recently been discussed in [27]. 2 More precisely usually variants of the equivalent photon approximation [32] are used [7,33]. In particular in the region of low ALP masses this approximation is problematic and receives sizable corrections that lower the cross section [27]. This further increases the relative importance of the production mechanism that we will discuss in the present paper.  [7] updated with the PrimEx recast [36].
Although the importance of the inclusion of axion production from secondary pions, was, for example, already pointed out in [34], it is still somewhat under-appreciated. Indeed all the proton beam dump constraints shown in the overview Fig. 1 are calculated using only the PFP production mode [7,33]. To our knowledge, only 3 a study put forward in [37] which determined prospects for SeaQuest (after its proposed ECAL upgrade) takes into account an estimate of the ALP yield stemming from the Primakoff-conversion of photons from π 0 decays in the dump 4 . However, some simplifying assumptions are made and no full Monte Carlo is set-up.
With our work, we want to close this gap in the literature and give improved estimates for the ALP production from meson-decay photons. The sensitivity improvement with respect to the case when only the PFP production is included will be discussed.
Our study is particularly timely since some of the experiments that can impact the parameter space have started taking data, notably NA62 [38], or are close to data-taking.
In practice calculating the photon flux inside the beam dump is far from trivial due to the non-perturbative nature of meson production. We therefore first carefully compare the yields for π 0 and other meson and their related angular distributions from PYTHIA [39] simulations to data from past experiments and then use it to determine the photon flux inside the dump.
While here we are interested mostly in ALP production, the π 0 and other meson spectra in the dump are not only of relevance to the production of ALPs but can also be the source of Dark Photons and other exotic particles. Thus our work in comparing the yield from PYTHIA simulations to data is of more general interest.
Our paper is structured as follows. A first direct comparison of the PYTHIA simulation output with the experimental data is performed in Section 2: proton-proton and protonberyllium interactions are separately analyzed. In Section 3 we review and discuss the computation of the ALP yield through production from the meson-decay photons. Finally, in Section 4 we re-evaluate existing experimental constraints and make estimates for future sensitivities taking the additional production mechanism into account. We discuss the conclusions in Section 5.

Neutral meson yields in proton beam dumps
The simulation of the production rates of secondary mesons is a challenging task, since the formation of mesons is complicated due to non-perturbative physics. In this section we therefore validate our PYTHIA simulations with experimental data. While no measurements of inclusive neutral meson production are available at exactly the desired energies and target materials employed in the experiments considered in Section 4, we nevertheless have data covering the energy range from 60 GeV to 450 GeV and different target materials, in particular hydrogen and beryllium. Furthermore we can compare to the production of different types of mesons. Putting this together allows us to have at least some confidence in the employed meson spectra. The total cross section for inclusive π 0 production agrees with the data within an uncertainty of 20% in the entire beam momentum range of interest.
We also indicate kinematic regions where the results are more uncertain. Including and neglecting the contributions from these regions we provide an estimate of the uncertainty of the limits and sensitivities in Section 4.

400 GeV proton beam on a hydrogen target
Measurements of secondary meson production from a 400 GeV proton beam dumped onto a hydrogen target have been performed at the beginning of the 1990's by the NA27 experiment operating at the LExan Liquid hydrogen Bubble Chamber (LEBC) with the European Hybrid Spectrometer (EHS). Results for π 0 and η production from the LEBC-EHS [41] allow a direct comparison of the proton-proton interaction expectation from the PYTHIA simulation program [39] with the experimental results. The simulation includes elastic, inelastic non diffractive, and single-, double-diffractive processes. Parton densities for protons are defined using the CTEQ 5L set [40], a widely-used leading-order QCD parametrization with α s (M Z ) = 0.127.
The measurements from LEBC-EHS report that π 0 (η) are produced with an average multiplicity of 3.87 ± 0.12 (0.30 ± 0.02) per incident proton. A total of 51.2 ± 3.1 % of the produced π 0 s stem from the decay of secondary particles (mostly mesons). These figures can be compared with the output of the PYTHIA simulation. From it, the total π 0 (η) production multiplicity is 4.248 ± 0.007 (0.489 ± 0.002). A total of 48.1 ± 0.7 %  of the produced π 0 s stem from the decay of ρ ± , ρ 0 , ω, or η, in good agreement with the experimental data. For the total cross section, the SoftQCD set of PYTHIA version 8.2 yields 39.9 mb summing up single-and double-diffractive, non-diffractive, and elastic processes, in good agreement with the experimental data which slightly exceeds 40 mb [42]. In Figs. 2 and 3, the distributions measured at LEBC-EHS for the squared transverse momentum (P 2 T ), rapidity (Y ), and Feynman X F variable have been compared to the PYTHIA output, after applying to the Monte Carlo the experimental selection criteria: for π 0 (η), the condition is X F > 0.006 (0.021), where the Feynman variable is computed as s, Z represents the beam axis direction, and P * is evaluated in the center of mass frame. To obtain the differential cross sections, the MC output is scaled according to the total proton-proton cross section, 39.14 mb for a 400 GeV proton beam. Four regions of X F are defined: a central region, for X F < 0.025, where the MC overestimates the data by a factor less than two; an intermediate region, for 0.025 < X F < 0.1, where data and MC agree; a fragmentation region, for 0.1 < X F < 0.3 where MC underestimates the data by a factor less than three, and a forward contribution, for X F > 0.3, associated to the inelastic diffraction mechanism, where the MC largely underestimates the data. Inserting in the simulation the double-pomeron exchange with the Minimum-Bias Rockefeller model [43] is seen to slightly improve the data-MC comparison for the central region, while not affecting the other regions. Including the double-pomeron exchange increases the total proton-proton cross section by 0.47 mb for a proton beam energy of 400 GeV.
In the following, the X F domain has been divided into the eight bins defined in Tab. 1. In general, the probability that photons from π 0 decays produce ALPs in the experimental acceptance increases with X F . As discussed in Section 3, we use our knowledge of the quality of our simulated meson spectra by making the following estimate for the uncertainty. As a baseline we take into account the full simulation results including all X F bins (0-7). This can then be compared to a conservative estimate that only includes bin 5, where the agreement between simulation and data is very good, and bin 6 where the simulation underestimates the data by a moderate factor of up to 3.

400 GeV proton beam on a beryllium target
We can now take the next step and allow for different target materials. First we note our expectation that 400 GeV protons interacting with a fixed proton target do not differ significantly from those of a neutron target, as far as the π 0 or η production is concerned. This is confirmed by the PYTHIA Monte Carlo both for the scattering distributions and cross sections. Following from this we obtain the cross sections for larger target nuclei by Figure 2. Measurements from scattering of 400 GeV protons onto a hydrogen target for π 0 (red) or η mesons (blue), compared with the expectation from PYTHIA 8.2: differential cross sections as a function of the squared transverse momentum (left panel) and rapidity (right panel). Symbols refer to results from the NA27 experiment at the LEBC-EHC. The histograms are the output of a PYTHIA MC simulation.  Figure 3. Measurements from scattering of 400 GeV protons onto a Hydrogen target for π 0 (red) or η mesons (blue), compared with the expectation from PYTHIA 8.2: differential cross section as a function of the Feynman X F variable. Symbols refer to results from the NA27 experiment at the LEBC-EHC. The histograms are the output of a PYTHIA MC simulation.
an appropriate scaling with the geometric cross section ∼ A 2/3 , which in this subsection we compare to the data.
Measurements of meson yields from a 400 GeV proton beam on beryllium targets have been performed in [44], and data was taken for four values of the secondary particle momenta (60, 120 and 300 GeV) and two values of transverse-momentum (0 and 500 MeV) at different target lengths. To complement these measurements at a lower momentum range of secondary particles and in view of evaluation of neutrino fluxes for NOMAD and CHO-RUS, the NA56/SPY experiment [45] published yields in the range of secondary momentum 2) with the same acceptance cuts applied by the experiment. Symbols represent NA56 data [45] (black dots for π + , red squares for π − ), black and red histograms represent the PYTHIA output for a 450 GeV proton beam. The dashed blue histogram indicates the expected inclusive invariant cross section for π 0 emission. The lower black and red histograms from PYTHIA refer to K + (black) and K − (red) and can be compared to the NA56 results for K + (black up-pointing triangles) and K − (red down-pointing triangles). Righthand-side plot: Inclusive invariant cross sections for protons (antiprotons) from NA56 represented by black circles (red squares), to be compared to black (red) histograms from PYTHIA. from 7 to 135 GeV with a proton beam of 450 GeV. In order to make these experimental data useful for further applications, a very useful parametrization was developed in [46], and sometimes is referred to as 'BMPT' 5 . As this parametrization was developed by the extrapolation of 400 GeV and 450 GeV data, it is suited to be employed for the use-case of NA62 and SHiP [6,56], while care has to be taken, when extrapolating to the NuCal [57], and SeaQuest beam energies of 70 GeV and 120 GeV, respectively. If the contribution of heavier meson and resonance decays is removed, the π 0 inclusive cross section is expected to be approximately equal to the average of the π + and π − inclusive cross sections. To validate our estimates for π 0 based on PYTHIA, we directly compare the MC output to the inclusive invariant cross section E d 3 σ/dp 3 obtained in [45] for π + and π − . For completeness, we also consider the emission of K + and K − , protons and anti-protons. Figure 4 shows the results of our comparison, which is quite good. To arrive at this comparison we have 1. accounted for the angular acceptance of NA56, corresponding to emission angles below approximately 0.7 mrad; 2. accounted for the NA56 target-efficiency factor: 1 − exp(L/λ Be ) (where L is the target length and λ BE is the appropriate proton interaction length) 3. scaled the MC results for the proton-beryllium total cross section, a factor of A 2/3 higher than that for proton-proton scattering.
Above, we have assumed a target length L =100 mm and a proton interaction length λ Be =423 mm. In Figure 4, left-hand side, we observe a very satisfying agreement for all available data from NA56 except a slight underestimation of the yield of π + at large momenta. As NA56 data might include tertiary production up to a certain extent, an underestimation with respect to the simulation output might be expected. On the right-hand-side of Figure 4, for completeness we also show the proton and anti-proton inclusive invariant cross sections measured at NA56 compared to the expected output from PYTHIA. 6

Proton beam energies below 400 GeV
To compare the inclusive production of light mesons with available literature, PYTHIA simulations of proton-proton interaction have been produced for proton beam energies of 70, 120, and 250 GeV: the first two values correspond to the beam energy of the NuCal and SeaQuest experiments, while the third value corresponds to the NA22 experiment [48], providing the most complete experimental data available below 400 GeV. The total cross section of π 0 production and the average number of emitted π 0 mesons are shown in Fig. 5 as a function of the beam energy: a general agreement (within 20% relative uncertainty) is observed between data [49] and MC. Data and MC differential cross sections have been compared at 250 GeV, as a function of the squared transverse momentum P 2 T and of the Feynman X F variable (Fig. 6). Data and MC transverse momentum differential cross sections are seen to agree (even if the data range is limited), while MC underestimates the data for X F > 0.05. Again, this will lead   to conservative estimates of the expected ALP yields in the following. The correlation between average transverse momentum and X F observed in data is quite well reproduced by the MC, as shown in Fig. 7.

Production of ALPs from meson decay photons
Having validated our PYTHIA spectra for neutral secondaries produced in proton interactions against data in the previous section 2, we are now ready to compute their impact on ALP production. We will also compare it with the contribution of ALPs from the photonfrom-proton-mode, which has been evaluated in Ref. [33]. As both processes have different initial and final states there is no interference and both contributions can be be added together. In the following we will therefore concentrate on the new contribution from the mesons.
Taking the meson distributions as an input we set up a Monte-Carlo simulation that proceeds along the following steps: 1. Decay of the neutral mesons produced in the dump into photons.
2. Compute the cross-section of ALP production from these photons in the target nucleus.
3. Mimic appropriate experimental acceptances and cuts.
4. Evaluate a sensitivity prospect at fixed number of incident protons for the situation of zero background.
Let us start with the first step. The decay length of the neutral pion is given by The decay length of the other neutral mesons is even smaller. Therefore, effectively the meson decay into photons is instantaneous, i.e. all mesons decay inside the target. The meson decay therefore yields a distribution of real photons depending on the energy E γ and the angle θ γ with respect to the beam axis. Due to the symmetry there is no dependence of this photon distribution with respect to rotations around the beam axis. The distribution is shown in Fig. 8 for proton beam energies of 400 GeV and 70 GeV.
In the second step the photons produce ALPs via the Primakoff conversion on the nucleus. We use laboratory-frame coordinates, where the nucleus is at rest. By a suitable rotation we can choose the photon to be moving in the z-direction. We then have, Note, however, that in this coordinate system the z-direction is slightly rotated compared to the beam axis, with the angle arising from the angle of the meson production and subsequent decay into the desired photon. In the Monte-Carlo we then rotate back to the system with the z-axis along the beam axis 7 . We use the approximate form of the cross section as in [33] 8 . In this simple coordinate system the cross section, approximated for small ϑ, reads, As in [33] the electromagnetic form-factor F em , is taken to be of the Helm form 9 (cf. [51]), To simplify the evaluation we set the form factor to zero for values q R 1 ≥ 4.49, i.e. above the first zero of the Helm form factor (as in [33]). Let us also emphasize again that in the present case we are dealing with a distribution of real, on-shell photons. This is in contrast to the photon from proton (PFP) mode where we have an effective parton distribution of virtual photons. Therefore, in this production mode we do not expect to be affected by the corrections to the equivalent photon approximation discussed in [27]. 8 We have checked by comparing to more complete expressions in [36,50] that the approximation to the cross section is excellent in the regime giving relevant contributions within the experimental acceptance. 9 At very low momentum transfer q 10 keV (depending on the target material) the electron shell also shields the charge, reducing the form factor. However, as discussed in [33] this region only gives a very small contribution to the signal which we neglect here.
The cross section is then determined by folding the (probability) distribution of the photons from the mesons with the cross section for an individual photon 10 As an example we show in Fig. 9 the resulting cross sections for several example energies and target materials relevant for the analysis in the next section (solid lines). This is then compared to the cross section for the PFP mode shown as dashed lines. For all considered energies 70 GeV, 120 GeV and 400 GeV, the PFP mode is sub-dominant. Also, for all energies, production from mesons is particularly favored since the spectrum is also considerably harder. The harder spectrum is useful for detecting ALPs with relatively large couplings since the higher γ-factor allows ALPs to decay outside dump and shielding regions.
For a better understanding of the behavior shown in Fig. 9 we note that the Primakoff cross section, Eq. (3.3), is peaked at small angles. More precisely, it is peaked at angles ϑ ∼ m 2 a /ω 2 . For sufficiently small masses and high energies this angle is very small. The ALP then has the same angle with respect to the beam axis as the incoming decay photon. Therefore, to a large degree the energy and angular dependence is dominated by that of the initial photon distribution shown in Fig. 8. However, this is only an approximate statement since the cross section (3.3) also has a significant tail towards larger angles where it is ultimately cut off by the form factor.
To get a better quantitative feeling for the involved scattering angles let us note that the radii, according to Eq. (3.5), are in the range R 1 ∼ (2 − 5) fm for the elements we are considering. With |q| ∼ E a ϑ for small masses, and the form factor setting in at |q|R 1 ∼ 1 we find that the tail is typically starting to be cut off at angles ϑ ∼ (0.002−0.004) (20 GeV/E a ). For energies 20 GeV, initial photon angles within a few milliradians of the desired ALP angle contribute. Stated differently, for energies E a 20 GeV and angles 10 milliradians, the ALP angle is dominated by the angle of the decay photon given in Fig. 8. Fig. 9 also demonstrates our procedure to estimate the uncertainty in our production rate. For 400 GeV we show two lines for the meson production. One (purple) in which we include all mesons in the full kinematic region and a more conservative estimate (black) where we include only the part of the meson spectrum that is most trustworthy. To be more precise we define two regions in X F . One covering all the bins, 0-7, in Tab. 1, whereas the other, more conservative one includes only the bins 5-6 where we find better agreement between Monte Carlo and data in our validation procedure. While this is very conservative, it does not fully appreciate the degree to which we underestimate the cross section. Indeed bin 6 gives a significant contribution to the cross section. From Fig. 3 we can see that in this region the simulated cross section underestimates the cross section by a factor which can easily be ∼ 2. Our conservative estimate therefore underestimates the cross section by perhaps few × 10%. Validated spectra improving in this region could therefore allow for a significantly increased sensitivity. 10 Let us note at this point that formally the PFP production and the production from decay photons is at different order in α. Meson production is a strong process and therefore essentially independent of α whereas "radiating" a photon from the proton requires an extra electromagnetic interaction. This is accounted for by a factor of α in the relevant photon distribution functions for the proton (cf. also [37]).    Table 1. These curves are shown at a primary proton energy of 400 GeV in purple (X F bins 0-7) and black (X F bins 5-6), 120 GeV (red, X F bins 5-6) and 70 GeV (brown, X F bins 5-6). The chosen target materials are copper at 400 GeV, and iron at 120 GeV and 70 GeV.
In the third step we then take into account all the relevant experimental cuts, • The ALP decay has to happen outside the target and in the decay volume.
Crucially this gives the exponential dependence on the decay length, where ALP is the decay length, d the distance from the target to the decay volume and l the length of the decay volume.
• The decay photons of the ALP have to reach the detector with a suitable minimum energy and other criteria are required for the photons to be detected.
We will describe the relevant details and approximations when we discuss the individual experiments. After all those cuts we obtain the fiducial cross section, σ f = fiducial cross section after all cuts. (3.6) Before we take the final step, let us also note that we have employed the simplification that both the meson production and the photon conversion into ALPs happen at the beginning of the target. This neglects the finite distance traveled before the proton interacts inside the target as well as the additional distance traveled by the photon before it is converted into an ALP. The typical distances are of the order of the proton and photon radiation length, respectively. As both are of the order of cm in the relevant energy range, this should be a minor effect compared to the target sized of the order of m for the experiments we will consider in Section 4. For experiments with targets in the cm range this would have to be taken into account.
In the final step we now need to compare the fiducial cross section for the detectable ALP production from a photon with the cross section for the photon to be absorbed in the target material. We then get for the total number of events, where N γ is the total number of photons produced in the meson decays. σ γN is the total cross section for photons to be absorbed in the target material and can be determined from the radiation length rad target , σ γ,target = m N ρ target rad target , (3.8) where m N is the mass of the nucleus. We stress that this is different from the PFP mode, where the relevant cross section would be the cross section for proton-nucleus interactions.
For the photon cross section we use values for the radiation length from [53], whereas for the proton cross section we employ σ pN 53 mb × A 0.77 from [54], where A is the mass number of the nucleus. For example in the case of copper the photon cross section for energies 2 GeV is about 6 barn for copper, 5.1 barn for iron and 12.7 barn for molybdenum. This is in contrast to the the proton cross section which is only 1.3 barn, 1.2 barn and 1.8 barn, respectively.
While the photon absorption cross section is significantly larger we find that in many cases this is more than compensated by a number of other factors. In particular we produce on average more than 1 meson per proton and each meson gives two photons. Importantly, as we have seen in Fig. 9, the spectrum of ALPs from these photons is also harder than the one from the PFP mode, which is advantageous for the detection in particular at larger ALP masses.

Updated sensitivities for fixed-target experiments
In the following, we will first update the exclusion contours for the past fixed target experiments CHARM [55] and NuCal with ALPs produced from the decay of π 0 . We then project sensitivities for the existing NA62 [38] and SeaQuest [37] set-ups. Finally, we give projections for the proposed SHiP [6,56] facility.

Past experiments: NuCal and CHARM update
For CHARM [55], we make use of the following parameters: The detector was located at a distance d = 480 m away from the proton dump, and the protons were dumped into a copper target. The detector (for CHARM this is identical to the decay volume) was l = 35 m in length and 3 × 3 m in transverse dimensions. The detector is off-set transversally by 5 m from the beam axis and this is accounted for in the MC. According to Ref. [55], CHARM was sensitive to events with a single electromagnetic shower in acceptance. CHARM quotes a signal acceptance of 51%, which we include in our estimate. The number of protons on target (POT) is 2.4 × 10 18 .
NuCal [57], made use of the U70 proton beam facility with a beam energy of 'only' 70 GeV. However, NuCal profits from a comparably small distance between target and detector of only d = 64 m and a detector length of 23 m. We adapt the analysis strategy of [58], and require a minimum ALP energy of E a > 10 GeV and at least one photon detected. In this way, the acceptance is approximately constant and equal to 70%. The detector has a radius of 1.3 m. In a dataset of N pot = 1.7 · 10 18 protons on an iron target, NuCal observed 1 event compared to a background expectation of 0.3 events. At 90% confidence level, we can therefore exclude any point in the parameter space predicting more than 3.6 events.
The changes in the limits from CHARM and Nucal when including ALPs produced from decayed π 0 s, can be understood from Figure 9.
For NuCal, where the beam energy is 70 GeV, the plot in the right-hand side of the figure illustrates that, including the π 0 yield appreciably changes the existing limits. This can be seen in Figure 10 (l.h.s.), where the brown dashed line shows the NuCal limits in PFP mode only, while the yellow region is the reach considering the added yield for NuCal from PFP and a conservative (bins 5-6) X F range from decayed mesons. However, the position of the upper part of the exclusion contour (at large couplings) is mostly determined by the experiment's geometry. Thus, even the revised NuCal limit does not drastically alter the untested parameter space that should be probed by new experiments, discussed in the next section.
Also for CHARM, including the ALP production channel via π 0 s improves the sensitivity with respect to the CHARM reach from PFP only considerably 11 . Figure 10 (l.h.s.) shows a zoomed version of the existing limits with updated CHARM contour (solid red) compared to the previous curve (dashed magenta). However, the CHARM limits lie still within the limits of E137 and the new NuCal limit. Thus we summarize that the experimental future landscape, does not change significantly even after the CHARM and NuCal updates. However, as we will see is the next subsection, the inclusion of the meson contribution, drastically changes the prospect sensitivity for forthcoming searches.

Current and future set-ups
The NA62 experiment [38] has been built to achieve a precise measurement of the ultra-rare decay K + → π + νν. Besides its main goal, NA62 has a rich program to search for exotic particles, including long-lived particles that can be produced in the up-stream copper beam collimator, into which the primary SPS proton-beam is fully (in dump-mode) or partially (in standard, parasitic data-taking) dumped (see, e.g. [60] for more details).
To be sensitive to a fully neutral final state, NA62 has to be run in beam-dump-mode. In our MC, we model NA62 using the following parameters: The distance between the beam-defining collimator (used to dump the beam) and the start of the fiducial volume is d = 82 m and the vacuum decay region before the Liquid Krypton Calorimeter (LKr) is l = 135 m long. In addition, we require the following acceptance conditions: Both photons produced in the ALP decay need to be detected at a minimum mutual distance of 10 cm in the LKr. Moreover, these photons need to be 15 cm away from the LKr central hole and their combined energy needs to be above 3 GeV. The target material for NA62 is copper and we show NA62 prospects for two different choices of POT.
In Figure 10 (l.h.s.) we show the sensitivity for NA62 at 1.3×10 16 POT (corresponding to a one-day run), using (purple dashed) PFP only, using the full X F range (bins 0-7, blue dotted) and the central X F bins 5-6 (red), respectively. We also show the sensitivity for 1.3 × 10 16 POT PFP to facilitate comparison with previous results [33]. As shown, the contribution from the decay of π 0 's dominates the sensitivity prospect for a 400 GeV proton beam and constitutes one of the most relevant results of our study. Regarding the question, whether all or only a conservative number of X F bins should be chosen, we observe that the difference between the reach in both situations is small but visible at ALP masses of few tens of MeV and small couplings. This can be understood by looking at Figure 9. The central π 0 production corresponding to X F bin 4 contributes to low photon momenta and correspondingly low ALP momenta. As we can see in Fig. 10 this is more relevant for the sensitivity at low masses and couplings. The forward contribution to π 0 production corresponding to bin 7, is not relevant in increasing the cross section at high photon energies. Not including X F bins other than 5 and 6 has a minor impact on our projections and thus for all following computations we adopt this conservative condition.
In Figure 10 (r.h.s) we show again the prospect for NA62, this time however at 10 18 POT (corresponding to a few months of data taking) while summing up contribution from π 0 production and PFP. Comparing this to the prospects shown, e.g. in [15,60] underlines the importance of including the yield of ALPs produced by π 0 s in these estimates.
The sensitivity reach for SeaQuest to ALPs produced as a result of secondary π 0 decays was estimated previously in [37]. As outlined in [37], the sensitivity to a di-γ final state requires the installation of an ECAL, potentially adapted from the PHENIX detector at BNL. Our analysis improves the study put forward in [37] in several key regards. Firstly, as outlined in Sect. 2, we have validated the π 0 differential cross sections obtained in PYTHIA against experimental literature in a wide energy range. Secondly, as for all considered setups, we have implemented a full Monte Carlo of ALP production and decay according to the geometry laid out in [37]. Lastly, we consider also the sub-dominant PFP contribution for ALP production in our estimate.
We use the geometric setup described in [37]. Also, as in [37], we assume the need for 10 signal events to detect a signal beyond the background fluctuations. The calorimeter is placed between tracking stations 3 and 4 (at ∼ 18.5 m downstream of the target) We use a fiducial volume that has a length of 1 m, in between meters 7 and 8 of the experiment, and a geometric acceptance of 2 × 2 m in transverse directions. The target material of SeaQuest is iron and we assume the phase-I statistics of 1.44 × 10 18 POT. Similarly to other set-ups, we require both photons to be detected at a minimum energy of 1 GeV each, a total energy of at least 3 GeV and a minimum mutual distance of 10 cm to avoid shower overlap given the photon shower Moliere radius.
The resulting prospects for SeaQuest are shown as the brown curve on the r.h.s. of Figure 10. Compared to the estimates of [37], the expected sensitivity covers a somewhat larger area of parameter space. Finally, we model the prospects for detection of ALPs in the SHiP [6,56] calorimeter as follows 12 : The fiducial region is taken to be 45 m downstream of the production point. The calorimeter is positioned at 50 m after the beginning of the decay volume. We ask both photons to be in an acceptance area of 5 ×10 m 2 . Both photons should have a minimum energy of 1 GeV and a combined energy of 3 GeV and be at least 10 cm apart. The result can be seen in Fig 10 r.h.s as green curve. Target material is molybdenum and the POT are 10 20 . Note that the envisaged SHiP calorimeter [62] has the potential of reconstructing the photon direction, thereby allowing ALP mass reconstruction. Compared to the results shown in [15] the mass reach increases considerably, from ∼ 1 GeV to ∼ 1.5 GeV.

Conclusions
Proton beam dump experiments are a popular and versatile tool to explore the dark sector in the MeV to GeV range that may be connected to a number of open problems in particle physics, notably dark matter. Important examples of current and near future experiments are NA62 (in beam dump mode), SeaQuest, and SHiP.
New particles can be produced in primary interactions of the proton with the target material but also in the decay of secondary mesons. However, further important contributions to the production can arise from the interaction of secondary or even tertiary particles with the target material. While these production mechanisms have been noted and even occasionally used [35,37] they are still somewhat under-appreciated and many sensitivity calculations do not take them into account.
In this paper we have performed a detailed investigation of axion-like particle (ALP) production. More precisely we discussed the production of ALPs from the following process: Protons interact with the target nucleus and produce neutral mesons. The mesons (mostly π 0 ) decay into two γ which subsequently can interact with another target nucleus to produce an ALP via the Primakoff process. We show that this gives a significant contribution to the production of ALPs which is also kinematically well suited for detection in typical experimental setups. Indeed, for experiments with high beam energies such as NA62 or SHiP, this is the dominant contribution in the region of interest and significantly extends the mass reach, e.g. by a factor of ∼ 1.5 in the case of the SHiP experiment.
A crucial input for the calculation of the production with secondary or even higher order particles are the spectra of these particles inside the target. In particular for mesons theoretical predictions are challenging due to the non-perturbative nature of the meson production processes. We have therefore validated our simulation results from PYTHIA 8.2 against a variety of measurements, thereby giving an estimate of the reliability of the simulation results and the impact this has on the sensitivity calculation for the experiments. We find that the impact of the uncertainty is moderate despite a relatively large uncertainty of our generated meson spectra in some regions of phase space. Further, very desirable improvement could come from two directions. First of all, our simulations only include mesons produced from the primary proton beam but, also for meson production, secondary interactions may play a sizable role. Including these secondaries will be an important next step. Second, the discrepancy of data and MC should be clarified by an extended study of existing data or, if needed, new measurements.
All in all, interactions of secondary particles in the beam dump are a powerful additional production mode for new very weakly coupled particles. Further studies that go beyond the example presented in this work are needed and in preparation.