Low-mass dark sector searches with deuteron photodisintegration
Cornelis J. G. Mommers Marc Vanderhaeghen







cmommers@uni-mainz.de

Institut für Kernphysik and PRISMA^+ Cluster of Excellence, 
 Johannes Gutenberg-Universität, D-55099 Mainz, Germany





Recent years have seen much activity in searches for dark-sector messenger particles in the 10–100 MeV mass range, especially in view of a potential new light boson conjectured by the ATOMKI collaboration, X17.
Under the assumption that the messenger particle has definite parity and either zero or unit spin, quite stringent bounds already exist on its coupling to electrons and protons. 
Equally stringent bounds on the neutron coupling do not exist yet, but are nonetheless desirable.
We explore how measurements of deuteron photodisintegration with a quasi-free neutron can yield bounds on the neutron coupling, and compute projections for a potential measurement at the low-energy high-intensity electron scattering experiment MAGIX@MESA.
The projected bounds are found to be competitive for an axial-vector or pseudoscalar scenario, but not for a vector or scalar scenario.


Low-mass dark sector searches with deuteron photodisintegration
    Marc Vanderhaeghen
    April 29, 2024
===============================================================







§ INTRODUCTION
 

In recent years, the main paradigm for dark-matter searches has been shifting.  
With heavier candidate particles, such as WIMPs, either excluded or experimentally out of reach <cit.>, a vigorous effort is presently underway to search for light dark-sector particles in the MeV–GeV mass range <cit.>. 
In this context a new boson X is often posited to be a messenger particle connecting the known visible sector to a hidden dark sector <cit.>.

The mass of the mediator X determines its decay. It may either decay invisibly into particles of the hidden dark sector, or it may decay visibly into particles of the Standard Model. 
Herein we deal with visible decays of the mediator X to an electron-positron pair, where one typically searches for X by directing a high-intensity, 𝒪(100 MeV–few GeV) electron beam at a fixed gas-jet target <cit.>.
If present, a potential signal will appear as a narrow resonance on top of a large, but well-known QED background. 

Using these techniques robust exclusion limits on the coupling of X to an electron-positron pair have been established, with reactions such as e^- p → e^-p (X → e^+ e^-), or e^- Z → e^- Z (X → e^+ e^-), where Z denotes a heavy target nucleus <cit.>. 
Likewise, results from the NA48/2 experiment, using π^0 →γ(X → e^+ e^-), have put strict bounds on the coupling of X to the proton <cit.>. 
And, under the assumption that X is a pseudoscalar, results from the SINDRUM collaboration <cit.>, which measured π^+ → e^+ ν_e (X → e^+ e^-), limit the isovector nucleon coupling of X.

Unlike the effective electron and proton couplings of X, the effective neutron coupling is not as stringently tested in the region m_X ∼ 10–100 MeV/c^2. 
The reason is simple: in the lab, no free, high-density neutron target exists.
Any extractions of the neutron coupling have to be made by proxy, which is why the current strongest constraints come from neutron star phenomenology <cit.> or isotope shift spectroscopy <cit.>.
Unfortunately, these bounds fall short of covering the MeV-to-GeV mass range of interest. 

Nevertheless, for several reasons it is desirable to have better limits on the neutron coupling in this region as well. First, taken with the existing constraints on the proton coupling, constraints on the neutron coupling can be directly translated into limits on the quark couplings of X, i.e. g_u = (2 g_p - g_n) / 3 and g_d = (2 g_n - g_p) / 3.
In turn, these limits directly tell us what new physics models are permissible. As a specific example consider the dark photon which, through kinetic mixing with an ordinary photon, necessarily couples to the quarks in a manner proportional to their charges, g_u = 2 g_p / 3 and g_d = - g_p / 3, implying g_n = 0 <cit.>.
Any deviation from this value already rules out the dark photon scenario. 
Another more immediate reason is connected to the potential discovery of a narrow resonance around 17 MeV by the ATOMKI collaboration—dubbed X17—in the e^+ e^- decay spectrum of excited states of ^8Be, ^12C and ^4He <cit.>. 
Limits on the coupling of X17 to the nucleon, which can be derived individually from each of the observed nuclear decays, appear to be in tension with each other <cit.>. 
Better constraints on the neutron coupling would help resolve this apparent conflict.

The problem, which involves determining the X coupling to the neutron, is mirrored in the program to measure the neutron polarizabilities <cit.>. 
There, a successful strategy is to measure deuteron photodisintegration under specific kinematic conditions where the bound neutron behaves as though it were free <cit.>. 
One can then access the polarizabilities indirectly, using the `quasi-free' bound neutron <cit.>.
The purpose of this paper is to explore the use of a similar approach, a measurement of γ d → e^+ e^- pn in quasi-free neutron kinematics, to extract limits on the neutron coupling of X.

To make our analysis concrete, our work will be centered around the forthcoming experiment at MAGIX@MESA <cit.>, though it may be readily adapted to meet the constraints provided by other electron scattering facilities. 
MESA is an electron accelerator (presently under construction) that will be capable of producing a low-energy, high-intensity electron beam (up to 105 MeV in its energy-recovering mode).
With a gas-jet target, luminosities around 10^35 cm^-2s^-1 will be reached <cit.>.
The MAGIX experiment consists of two high-resolution electron-positron spectrometers with a relative momentum resolution of δ p / p < 10^-4 and an angular resolution of δθ < 3 mrad <cit.>. 
Accordingly, the expected electron-positron invariant-mass resolution is δ m_ee≈ 0.1 MeV/c^2 or better <cit.>.
Note that earlier dark-photon searches have already achieved a resolution of 0.2 MeV/c^2 <cit.>.

In this paper we will investigate which limits on the neutron coupling to a dark-sector messenger particle X can be obtained with the MAGIX@MESA setup using γ d → e^+ e^- pn. Our paper is structured as follows: in Sec. <ref> we start by specifying the kinematics of the γ d → e^+ e^- pn process and derive an expression for the differential cross section. 
In Sec. <ref> we subsequently describe the plane-wave impulse approximation, and outline the kinematic regime in which the scattering on a quasi-free neutron is valid. 
In Sec. <ref> we specify the QED background and X-boson signal, respectively.
In Sec. <ref> we derive projections for the reach at MAGIX@MESA, and in Sec. <ref> we apply our framework to the specific case of the conjectured X17.
We conclude with a brief outlook in Sec. <ref>. 



§ KINEMATICS AND CROSS SECTION OF Γ D → E^+ E^- PN


We begin our calculation with the kinematics of the γ d → e^+ e^- pn process. 
We work within the rest frame of the deuteron and set up our coordinate system such that the z-axis is along the direction of the incoming photon's three-momentum, and the neutron and photon three-momenta span the x-z plane.
We have,

    γ(q, λ)   d(p_d, M)    → e^+(p_+, s_+)   e^-(p_-, s_-) 
              p(p_p, s_p)   n(p_n, s_n),

where λ is the polarization of the incoming photon, s_± are the e^± helicities, and M, s_p and s_n are the deuteron, proton and neutron spin projection along the z-axis. 
In the lab frame, the four-momenta are given by

    p_d^μ   = (m_d,        0,   0,   0),   
        q^μ  = (E_γ,   0,   0,   E_γ) ,
     
        p_i^μ   = (E_i,   |p_i| p_i),   for  i = ±, n, 
        
        p_p^μ   = (E_p,  p_p),

where

    p_i    = (sinθ_i cosϕ_i,  sinθ_i sinϕ_i,  cosθ_i),  for  i = ±, n, 
    p_p    = q - p_- - p_+ - p_n.

We denote the masses of the nucleons, leptons and deuteron as  m_N, m_e and m_d, respectively. 
It is also useful to introduce the four-momentum q' = p_+ + p_- and the squared invariant mass of the dilepton system, q'^2 = m_ee^2.

We choose our kinematic variables as E_γ, |p_±|, the polar angles θ_± and θ_n, and the azimuthal angles ϕ_± and ϕ_n. 
It follows that the magnitude of the neutron three-momentum is given by,

    |p_n| = 1/2a[ b ±√(b^2 - 4ac)],

where

    a    = (E_γ + m_d - q'^0)^2 - |q - q'|^2 cos^2θ_nγγ > 0  ,
    
        b    = (q + p_d - q')^2 |q - q'| cosθ_nγγ ,
    
        c    = m_N^2 (E_γ + m_d - q'^0)^2 - 1/4(q + p_d - q')^4 > 0 ,

with

    | q - q' | cosθ_nγγ = p_n · ( q - q' ).

For our kinematic regime of interest the positive solution of Eq. (<ref>) corresponds to physical values of the neutron three-momentum. 

Lastly, the differential cross section is given by

    dσ/dΠ = K ⟨ | ℳ |^2 ⟩,

with

    K    = 1/ 64 ( 2π )^8 m_d E_γ | p_+ |^2 | p_- |^2 /E_+ E_-
           × | p_n |^2 /| p_n | ( m_d + E_γ - q'^0 ) - E_n | q - q' | cosθ_nγγ,

where dΠ is shorthand for d |p_+| d | p_- | dΩ_-  dΩ_+  dΩ_n. Furthermore, in Eq. (<ref>) |ℳ|^2 is the squared Feynman amplitude and ⟨ · ⟩ denotes the average over initial helicities or spin projections and sum over final helicities or spin projections.



§ THE PLANE-WAVE IMPULSE APPROXIMATION




As a first step, to calculate the matrix element ℳ we work within the plane-wave impulse approximation (PWIA) <cit.>.
In this approximation, shown in Fig. <ref>, we can separate the process γ d → e^+ e^- pn into two parts: a part where the proton is a spectator and a part where the neutron is a spectator. 
By doing so we disregard meson exchange currents and final state interactions. 
However, in our kinematic regime of interest the meson exchange currents are estimated to give corrections of approximately 5%, meaning they can be safely neglected.
Likewise, for a first approximation, the final state interactions can be omitted (see Fig. 5 and its discussion in Ref. <cit.>).

As we are interested in the low-energy regime (E_γ∼ 100 MeV) where relativistic corrections are expected to be small, we use a non-relativistic framework to implement the PWIA. 
We start by inserting a complete set of two-particle states [We do not need to include a label for isospin, as a potential isospin factor of 1/√(2) is already included in the CB-Bonn wave function <cit.>.],


    ℳ_IA(γ d → e^+ e^- pn)    = ∑_s_1, s_2∫d^3 p/(2π)^3 ⟨ e^+ e^-,  p_p   s_p,  p_n   s_n  | ℳ_IA | γ,  12p_d + p  s_1,  12p_d - p  s_2 ⟩
               ×(2m_d)^1/2/(2E_p̃)^1/2 (2E_ñ)^1/2 _NR⟨p;   s_1   s_2  |  d(1, M) ⟩_NR,


where `NR' indicates the baryon states are normalized non-relativistically and no subscript refers to covariant normalization,

    |p,   s ⟩   = (2 E_p)^1/2|p,   s ⟩_NR, 
    ⟨p,   s  | p',   s' ⟩   = (2 E_p) (2π)^3 δ^(3)(p - p') δ_ss'.


In the second line of Eq. (<ref>) we can identify the components of the relative deuteron wave function in momentum space, Ψ̃^M_s_1 s_2(p) = _NR⟨𝐩;    s_1   s_2  |  d(1, M) ⟩_NR. 
They are given by

    Ψ̃^M_s_1 s_2(p)    = (2π)^3/2[ ψ̃_0(p) Y^0_0(p) ⟨12  s_1;  12  s_2  |  1   M ⟩
       - ψ̃_2(p) ∑_M_s Y^M - M_s_2(p) ⟨ 1   M_s;   2   M - M_s  |  1   M ⟩
         ×⟨12  s_1;  12  s_2  |  1   M_s ⟩] ,

where p = |𝐩|, ⟨ j_1   m_1;   j_2   m_2  |  j   m ⟩ are Clebsch-Gordan coefficients and Y^m_l are the spherical harmonics.
For our numerical estimates we use the CD-Bonn parametrization <cit.> for the s- and d-wave functions. Note that we use a different convention for the Fourier transform as compared to Ref. <cit.>, which gives rise to the different factors of 2π.
Also note that the relative sign between the s- and d-wave functions is absorbed into the CD-Bonn parametrization of Ref. <cit.>. That is, ψ̃_0(p) = ψ̃^Bonn_0(p) and ψ̃_2(p) = -ψ̃^Bonn_2(p). 
The wave function is normalized as

    1 = ∫d^3 p/(2π)^3 | Ψ̃^M (p) |^2 = ∫_0^∞dp   p^2 [ ψ̃_0^2(p) + ψ̃_2^2(p) ].


In this work we will restrict ourselves to kinematics inside the neutron quasi-free peak (NQFP) <cit.>. 
Qualitatively, this means we only consider configurations in which the bound proton is a spectator, i.e. has low momentum. 
As a consequence the incoming photon primarily scatters off the bound neutron ñ, and contributions to the cross section from processes involving the bound proton are negligible.   

A more quantitative definition of the NQFP may be derived using the relative and total momenta of the bound neutron and bound proton. 
In the lab frame the total momentum is just the deuteron momentum, p_d = 0, and the relative momentum is the momentum of the outgoing proton, p = p_p = - p_ñ.
Near the NQFP | p | → 0, meaning we approach the singularity of the deuteron wave function,

    1/m_d - E_p - E_ñ = 1/2 m_N - Δ - 2 E_p,

where Δ≈ 2.2 MeV is the deuteron binding energy.
Close to this singularity, ñ is `nearly on shell', and so we define the NQFP by <cit.>

    E_p ≲ m_N + Δ/2 |p_p| ≲√(m_N Δ)≈ 45.7  MeV/c.


Expanding Ψ̃^M in Eq. (<ref>) in the NQFP region gives

    ℳ_IA(γ d → e^+ e^- pn) ≈ℳ^n_IA(γ d → e^+ e^- pn),

with

    ℳ^n_IA(γ d → e^+ e^- pn) = (2 m_d)^1/2(E_p/E_ñ)^1/2
         ×∑_s_ñΨ̃^M_s_p s_ñ( p_p ) ℳ( γ ñ→ e^+   e^-   n ),

where the quasi-free neutron has momentum and spin projection -p_p and s_ñ, respectively. 
We have explicitly checked the size of the proton contribution to the cross section and found it to be negligible in the kinematic region corresponding with Eq. (<ref>), validating Eq. (<ref>). 
Note that in the NQFP kinematic region of interest here any off-shell effects are small, and can be safely neglected.



§ QED BACKGROUND PROCESSES AND SIGNAL PROCESSES
 

Having related the amplitude of γ d → e^+ e^- pn to the amplitude of γñ→ e^+ e^- n, we turn to the calculation of ℳ(γñ→ e^+ e^- n).
This calculation is composed of two steps, the calculation of the QED background amplitude and the subsequent computation of the signal amplitude. 



 §.§ QED background amplitudes


For photon energies around 100 MeV, i.e. below the pion production threshold, the QED background may be parameterized by a combination of the Bethe-Heitler process and Compton scattering process. 
The latter is described by the Born process, the π^0 t-channel exchange process and the non-Born contributions parameterized by the neutron electric (α_E) and magnetic (β_M) polarizabilities. 
These are all shown in Fig. <ref>, where we have omitted the crossed diagrams for the Bethe-Heitler and Born contributions.
All diagrams in Fig. <ref> are embedded in the quasi-free neutron blob of Fig. <ref>.
Technical details of the calculation of the QED background have been relegated to Appendix <ref> and a comparison of our calculated QED background and previous calculations is given in Appendix <ref>.

 



 §.§ Signal amplitudes

The X signal process is identical to Fig. <ref>b, with X replacing the virtual photon.
In principle, X may also contribute via the Bethe-Heitler process. 
However, to extract limits on the neutron coupling from γ d → e^+ e^- pn, we are only interested in processes where X is on resonance (that is, m_ee^2 = m_X^2). 
Consequently, any Bethe-Heitler contributions stemming from virtual X exchange are negligibly small and may be safely omitted.

We assume that X is either scalar, pseudoscalar, vector or axial vector (S, P, V or A, respectively) and that it couples to the fermions f, where f may be a neutron, proton or electron (f = n, p or e, respectively).
We can then generically parameterize its interaction using the effective interactions,

    ℒ_i = ∑_f g^i_f f̅Γ_i f · X,   for  i = S,  P,  V,  A,

where Γ_S = 1, Γ_P = iγ_5, Γ^μ_V = γ^μ and Γ^μ_A = γ^μγ_5. The dot in Eq. (<ref>) indicates all Lorentz indices, should they be present, are to be contracted. 
In this approach we remain agnostic to any models describing X;
constraints on the effective couplings g^i_f can always be matched to constraints on couplings from a UV-complete model by integrating out the corresponding degrees of freedom of the microscopic theory.

Using Eq. (<ref>) we find that the signal amplitude is given by

    ℳ_i = B_i/m^2_ee - m_X^2 + i m_X Γ_X,   for  i = S,  P,  V,  A,

where B_i = ε(q, λ) · H_i · C_i · L_i, with 

    C^(μν)_i    = -e g^i_e g^i_n [ δ_Si + δ_Pi -  g^μνδ_Vi
            - ( g^μν - q'^μ q'^ν / m_X^2 ) δ_Ai ], 
    
        L^(μ)_i    = u̅(p_-, s_-) Γ_i^(μ) v(p_+, s_+), 
    
        H^μ(ν)_i    = u̅(p_n, s_n) {Γ_i^(ν)(p_ñ + q + m_N)/(p_ñ + q)^2 - m_N^2Γ^μ_q 
            + Γ^μ_q (p_ñ - q' + m_N)/(p_ñ - q')^2 - m_N^2Γ_i^(ν)} u(p_ñ, s_ñ).

The photon nucleon vertex Γ^μ_q, is described in Eq. (<ref>), with e > 0.



 §.§ Signal averaged over a single invariant lepton-mass bin


As can be seen from Eq. (<ref>), the final cross section around the NQFP with m_ee = m_X will generally depend on the coupling of X to the neutron and electron, as well as the (narrow) width of X, Γ_X.
In an experiment we would not be able to set m_ee = m_X exactly—instead, we will collect all events falling in [m_ee - δ m_ee/2,   m_ee + δ m_ee /2 ].
If δ m_ee is sufficiently narrow, ∼𝒪(0.1 MeV/c^2), the QED background will remain nearly constant across this bin and, given that Γ_X ≪δ m_ee, this is also the only bin in which the signal will reside.

To calculate the signal averaged over a single bin we start by fixing δ m_ee = 0.1 MeV/c^2 (bin widths achievable by MAGIX@MESA).  
Moreover, on resonance we may safely neglect any interference terms between the QED background and the X-boson signal due to the narrowness of the resonance,

    dσ/dΠ≈( dσ/dΠ)_QED + ( dσ_i /dΠ)_X,   for  i = S, P, V, A,

where

    ( dσ_i /dΠ)_X := 1/δ m_ee∫^m_X + δ m_ee / 2_m_X - δ m_ee / 2dm_ee ( dσ_i /dΠ)_X,

is the signal averaged over a bin. 
We assume that m_X ≫Γ_X, so that in Eq. (<ref>) we may replace,

    |1/m^2_ee - m_X^2 + i m_X Γ_X|^2 →π/m_X Γ_Xδ(m^2_ee - m_X^2).

This gives

    (dσ_i /dΠ)_X 
        = 1/δ m_eeπ/2 m_X^2 Γ_X( K ⟨| B_i |^2 ⟩)|_m^2_ee  =   m_X^2,

where K is given in Eq. (<ref>) and B_i is given in Eq. (<ref>).

We can rewrite Γ_X using 

    ℬ(X → e^+ e^-) = Γ(X → e^+ e^-)/Γ_X.

A short calculation gives

    Γ_i(X → e^+ e^-)    = ( g_e^i )^2/π[ ( δ_Pi + δ_Vi) p_cm
            + ( δ_Si + δ_Ai) p^3_cm / m_X^2 ] β_i,

where β_S = 1, β_P = 1/4, β_V = (1 + 2m_e^2 / m_X^2)/6 and β_A = 2/3, with p_cm = √(m_X^2 - 4m_e^2)/2.

Using the above, dσ / dΠ no longer directly depends on Γ_X and g_e^i, and only depends on the electronic branching fraction, which appears as an overall scaling factor. Throughout this work we set ℬ(X → e^+ e^-) = 1. If the branching fraction is not equal to unity, any result can be easily rescaled by a simple multiplication.



§ NEUTRON-COUPLING REACH FOR MAGIX@MESA


Next, we present projections for the reach of the neutron couplings for an experiment at MAGIX@MESA.
We consider a signal significant if it is n_σ standard deviations over the QED background.
Then, the reach at the n_σ confidence level is given by <cit.>

    |g_n^reach| = [ σ_QED/σ_X|_g_n = 1n_σ/√(L ×σ_QED)]^1/2,

where σ_i are the integrated signal and background cross sections and we denote the integrated luminosity by L = ∫dt  ℒ, with ℒ the luminosity.
Equation (<ref>) implies an approximate scaling of the reach with the bin width like √(δ m_ee). 



The range over which we can integrate the signal and background cross sections is restricted by: a) the quasi-free neutron condition (<ref>), b) the requirement that m_ee falls into the correct bin (with width δ m_ee = 0.1 MeV/c^2), and c) geometric constraints from MAGIX@MESA <cit.>. 
The latter limit |θ_±| ≤ 165° and |θ_n| ≥ 5°.
Even with these constraints the accessible phase space is much larger than can be realistically measured in a single run. 
To that end, we consider several more restrictive detector settings.
Across all settings we fix

    |p_±|    ∈ [10.0,   100.0]  MeV/c, 
    ϕ_+    ∈ [150.0°,   210.0°],
    ϕ_-    ∈ [-30.0°,   30.0°], 
    θ_n    ∈ [-30.0°,   -5.0°] ∪ [5.0°,   30.0°],

and vary the other kinematic variables as outlined in Table <ref>. 

In all settings we only consider backward lepton kinematics because the signal-to-background ratio is considerably more favorable as compared to forward lepton kinematics. 
This is because in forward lepton kinematics the Bethe-Heitler part of the QED background increases, thereby significantly worsening the signal-to-background ratio.

We stress that the integration limits we provide are estimates.
Naturally, they are to be refined once the parameters for the MAGIX@MESA setup are established.
However, projections in Ref. <cit.> suggest our estimates are likely to be realistic.
Furthermore, should there be any significant changes then our analysis can easily be repeated. 
Finally, note that in practice the lepton spectrometers are restricted to in-plane kinematics only <cit.>. 
To cover a larger range of ϕ_± one can equivalently measure the neutron at an out-of-plane angle.  

4pt    

-4pt

We perform the integration using Monte Carlo integration with 10,000 integration points, at which results converge to within 1–2% of each other, which we deem to be sufficiently accurate for our purposes here.
We use rejection sampling to implement the NQFP condition (<ref>) and to ensure all events fall within the desired m_ee bin <cit.>.

The reach, assuming n_σ = 2, is shown in Fig. <ref>, for a scalar (a), pseudoscalar (b), vector (c) and axial-vector (d) X scenario.
To calculate the reach we have taken the projected luminosity of MAGIX@MESA's gas-jet target, ℒ = 10^35 cm^-2 s^-1 <cit.>, together with a beam time of 2000 h, giving an integrated luminosity L = 7.2 × 10^8 nb^-1.
Of course, results in Fig. <ref> can easily be rescaled using Eq. (<ref>) to account for different luminosities. 

It is evident that a higher E_γ correlates with an improved reach, and the best (lowest) reach is obtained in the 20–85 MeV/c^2 range with E_γ = 105 MeV.
To access masses ≲ 20 MeV/c^2 at E_γ = 105 MeV requires measurements in which the leptonic polar angles θ_± are further backwards than is geometrically feasible for the MAGIX@MESA setup.
To access the lower mass region, one instead has to lower E_γ.
As can be seen by comparing settings 1a–c, this unfortunately worsens the reach for the lower masses. 

To further explore the relation between the reach and our kinematic variables, in Fig. <ref> we show several projections of the phase space (sampled 10,000 times) assuming a vector-like scenario for setting 1c, with δ m_ee = 0.1 MeV/c^2 and m_ee = 25, 40 and 55 MeV/c^2. In Fig. <ref> we have, at each kinematic point, taken the inverse of Eq. (<ref>) using the differential cross section instead of the integrated cross section, yielding |reach|^-1. In addition, we have normalized each point to the maximum value of |reach|^-1 out of the 10,000 sampled points, denoted by max |reach|^-1. 
We have checked, and broadly our conclusions also hold for other parity assignments, bin sizes and settings.



In Fig. <ref>a one sees that a smaller m_ee corresponds to more backward leptonic polar angles θ_±. 
This is because m_ee is proportional to the relative angle between three-momentum of the lepton pair, meaning if m_ee is small, this relative angle has to be small as well. 
It is also evident that backward lepton kinematics give a better reach than forward kinematics.
The reason for this is, as previously mentioned, that in the forward regime the Bethe-Heitler contribution is enhanced, thus washing away the signal.

In Fig. <ref>b we have plotted the lepton momenta |𝐩_±|. Here one sees that higher momenta generally correlate to a better reach. In passing, let us also mention that projections in ϕ_± (not shown) are fairly evenly distributed.

Finally, in Fig. <ref>c we have plotted the proton momentum |p_p| against the neutron angle θ_n.
For all m_ee a clear pattern is visible: a lower protonic three-momentum corresponds to a better reach, which is not surprising. 
From Eq. (<ref>) we see that the reach is determined by the signal-to-background ratio as well as the absolute size of the integrated QED cross section.
The absolute value of the proton three-momenta enters the cross section in Ψ̃^M(p_p).
The most prominent contribution there, the deuteron s-wave function, is peaked at small momenta values. 
Thus, one obtains a larger cross section for small |p_p|, which is exactly reflected in Fig. <ref>.

As an aside, we remark that a cut on small |p_p| would not be advantageous.
We have checked that the gain in signal-to-background (the first fraction in Eq. (<ref>)) is smaller than the loss in phase space (the second fraction in Eq. (<ref>)) after placing the cut.
Therefore, the net result is a worse reach.
We discuss this further in Appendix <ref>.

In Fig. <ref> we have plotted the differential cross section of the background and signal as a function of θ_n. Here we assumed an axial-vector X scenario with m_X = 17 MeV/c^2, though similar conclusions hold for other parity assignments and masses as well. We fixed E_γ = 105 MeV, δ m_ee = 0.1 MeV/c^2 and the other kinematic variables to,

    |p_+|    = 65.7  MeV/c^2,   θ_+ = 165.0 °,   ϕ_+ = -180.0 °, 
    
        |p_-|    = 20.1  MeV/c^2,   θ_- = 165.0 °,   ϕ_- = 0.0 °, 
    ϕ_n    = 0.0 °.

The QED background is shown in red. In the above kinematics m_ee = m_X, leading to the enhanced signal shown in blue. For the blue curve, g_n^A = 7 × 10^-5, while the light blue bands correspond to the range g_n^A = (5–9) × 10^-5. The angular dependence shown in Fig. <ref> is typical for kinematics in the NQFP. The maximum of the peaked cross section corresponds approximately to the angle θ_n at which |p_p| reaches its minimum, assuming all other kinematic variables are kept fixed.



Lastly, in Fig. <ref> we have plotted the integrated background and signal cross section for an axial-vector X scenario as a function of m_ee, with m_X = 30 MeV/c^2 and δ m_ee = 0.1 MeV/c^2 using detector setting 1c.
Analogous results hold for different values of m_X, the bin width, and different parity assignments.
The signal height is proportional to |g^A_n|^2;
for illustrative purposes only we have chosen g^A_n = 10^-5. 

In an experiment, the strongly-enhanced signal would clearly be visible inside a single bin. 

We may compare the obtained reach with some available limits in the literature. As can be seen in Fig. <ref>a, for the scalar scenario limits from Ref. <cit.>,

    | g_n^S |× (m_X  MeV/c^2)^-2≲ 2.1 × 10^-5,

already exclude the accessible parameter space.

Similarly, for a vector-like X we have constraints from ^208Pb-n scattering <cit.>, limiting 

    | g_n^V |×|126/208 g_n^V + 82/208 g_p^V |× (m_X  MeV/c^2)^-4≲ 4.3 × 10^-10.

The intensive search program for the dark photon has yielded strong bounds on the proton coupling <cit.>, which we may reinterpret as bounds on a more general vector X <cit.>,

    |g^V_p| ≲ 2.5 × 10^-4.

Combining both constraints implies

    |g_n^V| × (m_X  MeV/c^2)^-2≲ 4.3 × 10^-6,

which also already excludes the accessible parameter space, as can be seen in Fig. <ref>c.

For a pseudoscalar X current limits on the nucleon coupling are given in terms of the isovector coupling <cit.>,

    | g^P_p - g^P_n | ≲ 1.2 × 10^-3.

Such an inequality may be easily satisfied if the proton and neutron coupling are of a similar order, so the purely-neutron exclusion limits one gets from γ d → e^+ e^- pn become relevant again. 

A similar situation arises for an axial-vector X, where, out of all parity assignments, the couplings are least constrained. The strictest bounds comes from the KTeV anomaly <cit.>,

    1.3 × 10^-10≤( g^A_p - g^A_n ) g^A_e/(m_X  MeV/c^2)^2≤ 5.2 × 10^-10,

where the upper and lower bounds are given in units of (MeV/c^2 )^-2.
The electron coupling can be estimated to be 

    g^A_e = ±(1.52 ± 0.31) × 10^-4,

which gives two possible bands in which the neutron coupling may fall (see Fig. <ref> in Sec. <ref> and Fig. 5 of Refs. <cit.>).
Moreover, just as the pseudoscalar coupling, one can fine-tune the proton and neutron couplings to largely evade the bound on their difference. 
Therefore, bounds solely on the neutron coupling would be a very welcome complement to the constraints given above.

On the other hand, one should also consider recent results from the isotope-shift analysis of Xenon <cit.> (and similar systems <cit.>).
At low energies both a pseudoscalar and an axial-vector interaction are described by a Yukawa potential in position space <cit.>. When included in isotope shift analyses, one obtains bounds on the product of the electron and neutron couplings. These bounds typically go up to 1 MeV/c^2 and are of 𝒪(10^-11). If we extrapolate them to our energy range, we get an order-of-magnitude estimate |g_e g_n| ≲𝒪(10^-9). Current bounds on the electron coupling are 𝒪(10^-4) <cit.>, which would mean the isotope-shift analyses would give limits on the neutron coupling of 𝒪(10^-5), much stricter than the ones we get in Figs. <ref>b and d. Of course, it remains an open question whether atomic methods can be improved enough to access the 10–100 MeV/c^2 mass range any time soon.

Next, we discuss how one may further improve the bounds in Fig. <ref>.
A possible approach would be to decrease δ m_ee. 
However, given that halving δ m_ee to 0.05 MeV/c^2 only gives an improvement of the reach with a factor ∼ 1/√(2), significant improvements from this area are unlikely.
As previously stated and demonstrated in Appendix <ref>, the reach can also be enhanced by increasing the measured phase space.
Given that the ranges given in Table <ref> already represent (roughly) the largest ranges that would be realistically possible for MAGIX@MESA to measure, the only remaining variable is the integrated luminosity.
As we already consider a fairly long experiment given our choice of 2000 hours, we turn to the luminosity.
It is unlikely that ℒ of the gas-jet target with a deuteron gas can be increased enough at MAGIX@MESA to make up the difference, as the luminosity enters as ℒ^-1/4 into the reach.
However, the gas-jet target does allow for different gases with more neutrons—for example Xenon, with around 80 neutrons. 
One naively expects the luminosity (and by extension the reach) to scale with the number of available neutrons.
To that end, interesting further research would be to investigate the coupling limits around the NQFP of a more neutron-rich system. 

In sum, deuteron photodisintegration gives competitive limits on the neutron coupling for a pseudoscalar or axial vector X, though it is likely that in the future methods such as isotope-shift spectroscopy will yield improved limits.
That being said, whereas other methods give limits on the product of g_n with another coupling, an extraction using γ d → e^+ e^- pn gives limits directly on g_n. Moreover, it is expected that by replacing the deuteron with a more neutron-rich system the reach can be improved further.



§ FURTHER RESULTS RELATED TO THE X17 ANOMALY



As mentioned in the introduction, recent results of the ATOMKI collaboration have garnered significant theoretical and experimental interest; in a series of experiments <cit.> the collaboration claims to have found evidence of a new, light boson dubbed X17.

The ATOMKI collaboration looked at internal pair creation in decays of excited ^8Be, ^4He and, recently, ^12C nuclei.
In all three cases, an anomalous bump was found in the distribution of the emitted electron-positron pair's relative angle, with a statistical significance consistently exceeding 6 σ (see Ref. <cit.> for a review).
Working with members of the ATOMKI collaboration, researchers at VNU University of Science claim to have independently verified the original ^8Be measurements <cit.>. 

In the Standard Model, nuclear transitions in which an e^+ e^- pair is emitted are mediated by electromagnetic interactions and are well understood. 
They are sensitive to new physics appearing at the MeV scale, and thus the ATOMKI collaboration attributes their anomaly to the as-of-yet unseen X17, with a reported averaged mass around 17.02(10) MeV/c^2 <cit.>. 
Assuming definite parity, the beryllium results indicate X17 can be a pseudoscalar, vector or axial-vector particle <cit.>, while the carbon results point to a scalar, vector or axial-vector particle <cit.>. 
Theoretically, based on the ^8Be decays, models for X17 have been developed for the pseudoscalar, vector and axial-vector cases <cit.> that attempt to explain the ATOMKI anomalies while conforming to existing exclusion bounds. 
In particular, according to the vector model put forward by Feng et al. <cit.> X17 additionally must be protophobic (couple weakly to protons) to meet existing bounds from the NA48/2 experiment <cit.>. 
Experimentally, a global effort is underway to scrutinize the results of the ATOMKI anomaly, with new experiments at facilities such as CCPAC <cit.>, MEG II at the PSI <cit.>, JLAB <cit.>, New JEDI <cit.>, among many others <cit.>.

We investigate how the neutron-coupling limits from deuteron photodisintegration discussed above can be of added value in the current search for X17.
To estimate the coupling of X17 to the nucleon we compare models by Alves and Weiner <cit.> for the pseudoscalar case, by Feng et al. <cit.> for the vector case and by Kozaczuk et al. <cit.> for the axial-vector case. 
To match their notation with ours we have 

    g^(0)_XNN   = (g^P_p + g^P_n) / 2,    g^(1)_XNN = (g^P_p - g^P_n) / 2,
    ε_f    = g^V_f / e , 
    
        a_f    = g^A_f,

for f = p, n and e.



We can constrain X17's couplings by using the branching fractions of the ^8Be and ^12C decays reported by the ATOMKI collaboration <cit.>,

    Γ_X/Γ_γ|_^8Be(18.15)   = 6(1) × 10^-6, 
    Γ_X/Γ_γ|_^12C(17.23)   = 3.6(3) × 10^-6.

as well as from the branching fraction of ^4He decays. The latter are more complicated with furthers details given in Refs. <cit.>. Because the couplings of X17 have been extensively discussed in other works we shall only give a brief summary.
For more details we refer the reader to the relevant papers. 

The ^8Be(18.15) state, which is predominately isoscalar, is isospin mixed with the ^8Be(17.64) state, which is predominately isovector.
In our analysis we parameterize this isospin mixing with an isospin-mixing angle, θ_1^+, and an isospin-breaking parameter, κ <cit.>. 
Following Ref. <cit.>, we take θ_1^+ = 0.35(8) °, whence κ = 0.681 <cit.>. 

For a pseudoscalar X17 scenario, results from the SINDRUM collaboration <cit.> put a strong bound on the isovector coupling, 

    |g^(1)_XNN| ≤ 0.6 × 10^-3.

By following the procedure described in Ref. <cit.> we derive bounds on the isoscalar coupling from Eq. (<ref>). 
This gives

    |g^(0)_XNN| = (2.0 - 5.4) × 10^-3,    |g^(1)_XNN| = (0.0 - 0.6) × 10^-3.

Clearly, from Fig. <ref>b our projected limits from deuteron photodisintegration for a pseudoscalar X17 are not competitive.

For a vector X17 scenario the constraint provided by the NA48/2 experiment <cit.> leads to the protophobia condition (the green band in Fig. <ref>a), 

    |ε_p| ≤ 1.2 × 10^-3.

We derive the remaining neutron coupling from the ^8Be data as outlined in Ref. <cit.>, and from the ^12C data using <cit.>

    Γ_X/Γ_γ|_^12C(17.23) = k/Δ E( 1 + m_X^2/2 Δ E^2) | ε_p - ε_n |^2,

where k = √(Δ E^2  - m_X^2). 


We find

    ^8Be:    |ε_n| = (1.1 - 1.7) × 10^-2 , 
    
        ^12C:    |ε_n| = (2.6 - 5.3) × 10^-3.

The projected limits from γ d → e^+ e^- pn, shown in Fig. <ref>, would have to be improved by at least a factor of two to test these couplings. This is also visible in Fig. <ref>a, where the γ d → e^+ e^- pn projections span the entire parameter space.

In passing, we remark that like Refs. <cit.>, we find some tension between the carbon, beryllium and helium results for a vector-like X17.
The neutron couplings only overlap when the uncertainty of the ^8Be results is increased to around 3σ (see the orange, red and pink bands in Fig. <ref>a). 
In view of this there is an ongoing discussion in the literature surrounding X17's couplings. 
There is a general consensus that, as it stands, one cannot explain the ^8Be, ^4He and ^12C anomalies simultaneously without introducing some form of tension (see Refs. <cit.> for more discussion thereon. From pion decays, the blue bands in Fig. <ref>a, one may obtain additional constraints which increase the tension even more). 
Our carbon-derived neutron couplings agree with Ref. <cit.>, which would indicate that one should critically re-examine the beryllium data; a similar finding is echoed in Ref. <cit.>. 
A possible resolution of the discrepancy would be if X17 simply is not a vector particle. 
If X17 exists, and if it were in truth an axial vector, then interpreting its axial couplings as vector couplings would, naturally, lead to erroneous results.

To derive couplings for the scenario of an axial-vector X17 we need its nuclear matrix elements.
For the beryllium transition we take the matrix elements as calculated in Ref. <cit.>,

    ⟨^8Be(g.s.) ||σ̂^(p)||^8Be(18.15) ⟩   = (-0.38 ± 2.19) × 10^-2, 
    ⟨^8Be(g.s.) ||σ̂^(n)||^8Be(18.15) ⟩   = (-10 ± 2.6) × 10^-2,

which may be used in <cit.>

    Γ_X |_^8Be(18.15) = k_X/18π[ 2 + (Δ E / m_X)^2 ] 
         ×|⟨^8Be(g.s.) || a_p σ̂^(p) + a_n σ̂^(n)||^8Be(18.15) ⟩|^2.

Because there are no bounds on the proton coupling or neutron coupling individually, the axial vector coupling is only loosely constrained. 
Even taking into account the additional KTeV anomaly (<ref>) (the green bands in Fig. <ref>b) still leaves a large part of the parameter space accessible. 

Unlike the beryllium decay, the required axial carbon matrix element, ⟨ ^12C(g.s.) || D̂_3^σ ||  ^12C(17.23) ⟩, with D̂_3^σ the corresponding spin dipole operator, is not known <cit.>.
This relates back to the aforementioned tension.
We view the computation of ⟨D̂^σ_3 ⟩ as a pressing open problem.
Its resolution would immediately tell us whether one has a true discrepancy between the carbon and beryllium results—there would be no overlap between the carbon and beryllium measurements in the parameter space of the axial-vector couplings—or whether one has merely misinterpreted the data. Order-of-magnitude estimates suggest that the ^8Be and ^12C data could be harmonized using an axial vector X17 <cit.>.

Following the analyses in Refs. <cit.> we summarize the aforementioned constraints in Fig. <ref> for the two most promising X17 scenarios, vector and axial-vector, with projected limits from Fig. <ref> overlaid on top. Note that, unlike in Fig. <ref>, here the shaded regions indicate allowed regions of the parameter space. Our ^8Be limits differ slightly from Refs. <cit.> as we use a different mixing angle θ_1^+. We have checked that when using their mixing angle we recover the their results. The actual value of the mixing angle has not been pinned down exactly (for example Refs. <cit.> use different values, see also the discussion in Ref. <cit.>), meaning the available parameter space may be slightly larger than shown in Fig. <ref>. For both the vector and axial-vector scenarios, changing the mixing angle effectively rotates the ^8Be bands around the origin. Especially for the axial-vector scenario limits from γ d → e^+ e^- pn would constrain to what extent this rotation is permissible.  

In summary, if a calculation of the axial transition matrix element of ^12C(17.23) to the ground state remains outstanding, constraints from γ d → e^+ e^- pn would be subdominant for a vector X17 scenario, but for what is at present the most likely X17 scenario, an axial vector, helpful to limit the available coupling range.



§ SUMMARY AND OUTLOOK

Low-mass dark sector searches are a very active area of research, fueled further by possible sightings of a new 17 MeV/c^2 boson, X17.
Compared to the electron and proton couplings, the coupling of a new boson in the 10–100 MeV/c^2 mass range to the neutron is poorly constrained.
To that end, we studied the efficacy of the deuteron photodisintegration process γ d → e^+ e^- pn at MAGIX@MESA in limiting the effective coupling of a new scalar, pseudoscalar, vector or axial-vector boson in the 10–100 MeV/c^2 mass range, using kinematics corresponding to quasi-free production on a neutron.
It turns out that for a scalar or vector scenario other methods yield better constraints on the neutron coupling.
However, for a pseudoscalar and axial vector scenario current bounds only limit the difference between the proton and neutron coupling. 
In this case, deuteron photodisintegration has the potential to yield limits directly on the neutron coupling.
Specifically, for the axial vector scenario our projections from γ d → e^+ e^-pn yield competitive constraints for the neutron coupling in light of X17.

Furthermore, besides providing exclusion limits on the neutron, via the reaction γ d → e^+ e^- pn one can also study the neutron polarizabilities, as has been previously done using γ d →γ pn <cit.>. Such a study would be a worthwhile extension of the present work.



§ ACKNOWLEDGMENTS

The authors thank S. Schlimme, D. Barducci and C. Toni for helpful communications.
This work was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation), in part through the Research Unit [Photon-photon interactions in the Standard Model and beyond, Projektnummer 458854507 - FOR 5327], and in part through the Cluster of Excellence [Precision Physics, Fundamental Interactions, and Structure of Matter] (PRISMA^+ EXC 2118/1) within the German Excellence Strategy (Project ID 39083149).





§ CALCULATION OF THE QED BACKGROUND





 §.§ Conventions

Our computation of the QED background starts with the QED vertex following from the QED interaction Lagrangian.
In our convention,

    ℒ_QED = - J_μ A^μ = e ψ̅γ^μψ A_μ,

The photon nucleon vertex Γ^μ_q, is defined as

    ⟨ N(p')  |  J^μ(0)  |  N(p)  ⟩
       = -e u̅(p') Γ^μ_q u(p) 
       = -e u̅(p') [ F_1(q^2) γ^μ + i/2 m_N F_2(q^2) σ^μν q_ν] u(p),

where N = p, n and J^μ(x) is the electromagnetic current operator, with q = p' - p.
For the Dirac and Pauli form factors, F_1 and F_2 respectively, we use the parametrization given in Ref. <cit.>. 



 §.§ Bethe-Heitler amplitude

The tree-level diagram for the Bethe-Heitler process is shown in Fig. <ref>a.
The amplitude is given by

    ℳ_BH   = ie^3/(q - q')^2 ε_μ(q, λ) L^μν_BH H^BH_ν,

where ε^μ is the photon polarization vector, with

    L^μν_BH   = u̅(p_-,s_-) 
        {γ^μ( p_- - q + m_e /-2 p_- · q) γ^ν
            + γ^ν( q - p_+ + m_e /-2 p_+ · q)  γ^μ} v(p_+, s_+),
    
        H^ν_BH   =  u̅(p_n,s_n) Γ_q - q'^ν u(p_ñ,s_ñ) .




 §.§ Born amplitude

The tree-level Born amplitude is shown in Fig. <ref>b, and its amplitude is given by

    ℳ_Born = - ie^3/q'^2ε_μ(q, λ) L_ν^Born H^μν_Born,

where

    L^ν_Born   = u̅(p_-,s_-) γ^ν v(p_+, s_+),
    
        H^μν_Born   = u̅(p_n,s_n) {Γ^ν_-q'( p_ñ + q + m_N )/( p_ñ + q )^2 - m_N^2 Γ^μ_q 
            + Γ^μ_q ( p_ñ - q' + m_N )/( p_ñ - q' )^2 - m_N^2Γ^ν_-q'} u(p_ñ,s_ñ),

with Γ^μ_q given in Eq. (<ref>).



 §.§ Pion-pole amplitude

For the pion-pole contribution we need the pion-photon and pion-nucleon vertices. 
The pion-photon vertex follows from

    ℒ_π^0 γγ   = e^2/4 F_π^0 γγ F^μνF̃_μνπ^0,

where F̃_μν = 12ε_μναβ F^αβ and F_π^0 γγ(0) = 1/(4 π^2 f_π), with F^μν the photon field strength and f_π≈ 92.4 MeV the pion decay constant. 
Our convention is ε_0123 = +1.
The photon-nucleon vertex follows from

    ℒ_π N N = g_π N NN̅π^a τ_a i γ_5 N,

where τ^a are the Pauli matrices, N = (p, n) and g_π N N = g_A m_N / f_π, with the axial-vector coupling from nuclear β-decay g_A = 1.27.
The pion-pole amplitude is shown in Fig. <ref>c, and is given by

    ℳ_π^0 = - F_π^0 γγ g_π N Ni e^3/q'^2ε^μ(q,λ) L^ν_π^0 H^π^0_μν ,

where

    L^ν_π^0   = u̅(p_-, s_-) γ^ν v(p_+, s_+), 
    
        H_μν^π^0   =  i/(q-q')^2 - m_π^2ε_μναβ q^α q'^βu̅(p_n, s_n) γ_5 u(p_ñ, s_ñ).




 §.§ Non-Born neutron polarizability amplitude

Lastly, we account for nucleon polarizability contributions following the formalism in Ref. <cit.>.
In short, one may parameterize the unpolarized Compton process, N(k) γ(q) → N(k')γ(q'), often subtracting the Born and pion exchange contribution, with 18 independent tensors T^μν_i and their corresponding coefficients B_i(q^2, q'^2, q · q', q · P), where P = 12(k + k'). 
For unpolarized real Compton scattering only two tensors contribute (see Appendix A in Ref. <cit.>),

    M^μν_RCS, spin avg.   = B_1(0, 0, q · q', q · P) T_1^μν
            + B_2(0, 0, q · q', q · P) T_2^μν,

where

    T_1^μν   = - (q · q') g^μν + q'^μ q^ν,
    
        T_2^μν   = (2 m_N ν)^2 ( -g^μν + q'^μ q^ν/k · q') 
            - 4 (q · q') (P^μ - q · P/q · q' q'^μ) (P^ν - q · P/q · q' q^ν),

with ν = (q · P) / m_N. 
The low-energy expansions for B_1 and B_2 are given by

    B_1(0, 0, q · q', q · P)    = 1/αβ_M1,
    
         B_2(0, 0, q · q', q · P)    = - 1/α1/(2 m_N)^2 (α_E1 + β_M1),

where the values for the neutron polarizabilities are <cit.>

    β^n_M1   = (3.7 ± 1.2) × 10^-4 fm^3, 
    α_E1^n + β^n_M1   = (15.2 ± 0.5) × 10^-4 fm^3.


The polarizability amplitude is illustrated in Fig. <ref>d and given by

    ℳ_Pol. = -ie^3/q'^2ε^μ(q,λ) L^ν_Pol. H^Pol._μν ,

with

    L^ν_Pol.   = u̅(p_-, s_-) γ^ν v(p_+, s_+) 
    
        H_μν^Pol.   = u̅(p_n, s_n) M^RCS, spin avg._μν u(p_ñ, s_ñ).

where the blob in Fig. <ref> indicates the polarizability contribution.
Note that in our case √(q'^2)≪ 1 GeV, so using the real Compton amplitude in place of the single-virtual Compton amplitude is a very good approximation.



§ VERIFICATION OF THE QED BACKGROUND





In general, there is no calculation of the QED background of γ d → e^+ e^- pn that we can directly compare to. However, there has been a calculation of the QED background in γ d →γ pn <cit.>. Up to phase-space factors this is, in principle, identical to our calculation with the Bethe-Heitler contribution removed. Therefore, it can serve as a useful intermediate cross check. Of course, the methods used in Ref. <cit.> are not identical to our methods, so we cannot expect a completely identical result. Nevertheless, our results should agree to at least the same order of magnitude.

To that end, we have repeated the PWIA calculation of Ref. <cit.>. Using their notation, we have computed the in-plane differential cross section dσ / dΩ_γ'dΩ_n T_n as a function of T_n, where T_n is the kinetic energy of the outgoing neutron and Ω_γ' and Ω_n are the solid angles of the outgoing photon and neutron, respectively (see Ref. <cit.> for details of the calculation). This is shown in Fig. <ref> and should be compared to the dotted line in Fig. 5a of Ref. <cit.>. As can be seen, both calculations agree, giving a useful cross check of our QED background.



§ DEPENDENCE REACH OF NEUTRON COUPLING ON INTEGRATION RANGE




From Eq. <ref> it is clear that the reach of the coupling of X to the neutron depends on the size of the phase space we integrate over.
Given that the detector settings in Table <ref> are estimates, it is important to see to what degree changing the accessible phase space influences the reach. 
The ranges of |𝐩_±| are unlikely to change in the final MAGIX@MESA setup, and cutting the range of θ_n has little influence on the reach (see Sec. <ref>).
Therefore, we focus on the influence of θ_±, ϕ_± and E_γ.

The general trend is that increasing the available phase space improves the reach, albeit with diminishing returns.
This is illustrated in Fig. <ref>. 
Here, we take detector setting 2 (with m_X = 65 MeV/c^2 and δ m_X = 0.1 MeV/c^2) and either increase the integration range of ϕ_± or θ_± symmetrically, 

    ϕ_+    ∈ [180 ° - Δϕ_±,  180 ° + Δϕ_±] ,
    ϕ_-    ∈ [-Δϕ_±,  Δϕ_± ],

or

    θ_±∈ [120° - Δθ_±,  120° + Δθ_±],

respectively, keeping the other parameters fixed. We have assumed a vector-like scenario, however we have checked that the same conclusions hold for the other parity assignments.

In Fig. <ref>a we see that, as a function of E_γ, the reach has an optimum. The reason here is that shifting the energy roughly corresponds to translating the exclusion limits in Fig. <ref>—at a certain point one hits the kinematic boundaries.

In Fig. <ref>b we see that increasing the integration range to be more out of plane quickly increases the reach, before flattening out.
A similar trend is visible in Fig. <ref>c, when one increases the range of the polar angles.
So, if one integrates over a very small part of the phase space, it is worthwhile to try to access more.
If one already can integrate over a large area of phase space, increasing the phase space further is not as useful.



