Heavy Axion Opportunities at the DUNE Near Detector

As defined, the QCD axion has a gluonic coupling (for a recent review see DiLuzio:2020wdo )

through which it solves the Strong CP Problem dynamically. This can be seen explicitly by considering the QCD-generated axion potential Weinberg:1977ma ; DiVecchia:1980yfw ; diCortona:2015ldu

As the axion dynamically relaxes to its minima at $\langle a\rangle=-\bar{\theta}f_{a}$, it makes the effective $\bar{\theta}$ parameter vanish – solving the Strong CP Problem and explaining the smallness of the (as-yet unobserved) neutron electric dipole moment Abel:2020gbr . The above QCD-axion potential gives rise to the well-known relation diCortona:2015ldu ,

However, a variety of terrestrial, astrophysical and cosmological constraints (see e.g. Vysotsky:1978dc ; Raffelt:2006cw ; Cadamuro:2011fd ; Millea:2015qra ; Zyla:2020zbs ) require $f_{a}>10^{9}$ GeV for the QCD axion and therefore, the relation in Eq. (3) precludes observing an otherwise phenomenologically interesting, accelerator-observable parameter space where $m_{a}\sim 10$ MeV - $100$ GeV. Hence, it is interesting to ask whether there exist models solving the Strong CP Problem which can occupy this mass regime.

At the same time, from an ultraviolet (UV) perspective, various axion models often suffer from the so-called “Quality Problem” Kamionkowski:1992mf ; Barr:1992qq ; GHIGNA1992278 ; Holman:1992us . To see this, we recall that an axion can be realized as the Goldstone boson of a spontaneously broken Peccei-Quinn (PQ) $U(1)_{\rm PQ}$ symmetry. In the far UV, quantum gravitational effects are expected to break all global symmetries Kallosh:1995hi ; Banks:2010zn . Therefore, at low energies, the $U(1)_{\rm PQ}$ can, at best, survive as some accidental symmetry. To illustrate the severity of the Quality Problem, we consider a Planck-suppressed $U(1)_{\rm PQ}$-breaking operator,

In the above, the axion $a$ arises as the Goldstone mode of the (composite) field $\Phi\sim f_{a}e^{ia/f_{a}}$ at low energies, and the resulting minima-structure of the axion potential need not align with the QCD-generated potential. Therefore, unless the UV contribution in Eq. (4) is small compared to the QCD-generated potential in Eq. (2), the axion will relax to a minima dictated by Eq. (4) where generally $\langle a\rangle\neq-\bar{\theta}f_{a}$ and the axion solution to the Strong CP Problem will be spoiled. More concretely, given the constraint $f_{a}>10^{9}$ GeV for the minimal QCD axion, we see that unless we forbid all operators of the type in Eq. (4) up to $N=9$, the axion-solution to the Strong CP Problem no longer works. The question of why the UV theory should respect $U(1)_{\rm PQ}$ to such high quality is the Quality Problem.

In Ref. Hook:2019qoh a model addressing the Quality Problem was constructed, generalizing on previous work Rubakov:1997vp ; Berezhiani:2000gh ; Hook:2014cda ; Fukuda:2015ana ; Dimopoulos:2016lvn , in which the Strong CP Problem is solved through the presence of a $Z_{2}$-symmetric mirror sector (containing the primed fields)¹¹1For other approaches for addressing the Quality Problem, including some recent discussions, see e.g.  Kim:1984pt ; Randall:1992ut ; Choi:2003wr ; Fukuda:2017ylt ; Lillard:2018fdt ; Gavela:2018paw ; Cox:2019rro ; Alvey:2020nyh . The axion coupling is then given by,

This $Z_{2}$ symmetry is softly broken by the only relevant operator in the SM and the mirror sector, the Higgs masses. Consequently, the mirror Higgs VEV $\langle H^{\prime}\rangle$ can “naturally” be much larger than $\langle H\rangle$. In this case, the mirror quarks decouple at higher energies without impacting the RG running of the mirror QCD at lower energies. This results in the mirror confinement scale $\Lambda_{\rm QCD^{\prime}}$ being much larger than $\Lambda_{\rm QCD}$. Therefore, the QCD-axion potential receives a parametrically larger contribution from the mirror sector with Hook:2019qoh ,

At the same time, thanks to the $Z_{2}$ symmetry, this enhanced contribution still aligns with SM-QCD generated potential in Eq. (2). Thus the same axion solves the Strong CP Problem in both the sectors, and, since the axion potential is parametrically enhanced due to the presence of the mirror sector, it is less susceptible to the Quality Problem. Since the RG running of $\bar{\theta}$ happens at seven loops and the threshold effects happen at four loops Ellis:1978hq , the two $\bar{\theta}$ angles remain approximately equal in the IR despite the spontaneous $Z_{2}$-breaking. We now briefly summarize the theoretical constraints on our model in Fig. 2, while referring the reader to Ref. Hook:2019qoh for more detailed explanations as well as how to obtain a viable cosmology in this class of models.

Fig. 2 presents the theoretically-motivated region of parameter space for this heavy, high-quality axion model, as a function of the axion mass $m_{a}$ and the axion decay constant $f_{a}\equiv f_{G}/(4\pi^{2})$. Since our mechanism makes the axion only heavier, it can not populate the region labeled “QCD Axion” where it would be lighter than the QCD axion. In the region labeled “$f_{a}<\Lambda_{\rm QCD^{\prime}}$”, the axion EFT breaks down because $m_{a}>f_{a}$ in that region. In the region labeled “PQ Quality Problem,” the axion suffers the Quality Problem discussed in Eq. (4) due to operators with $N\geq 6$. Finally, in the region in the bottom right denoted “$H^{\prime}$ Quality Problem”, the mirror VEV $\langle H^{\prime}\rangle$ spoils the Strong CP solution via,

In particular, to avoid the Quality Problem from Eq. (7), we require $\langle H^{\prime}\rangle<10^{14}$ GeV, so there is only a maximal amount by which the axion can be made heavier in this scenario.²²2While a portal coupling $\lambda|H|^{2}|H^{\prime}|^{2}$ can be present, its primary effect would be to make the SM Higgs very heavy unless $\lambda$ is very small. We view this generically large contribution to the Higgs mass as another form of the hierarchy problem for the SM Higgs which we do not try to address in this work.

By inspection, Fig. 2 encourages us to focus on heavy axions in the keV-TeV mass range with $f_{G}$ between $10-10^{9}$ GeV. A natural question that emerges is how much of the open theoretical parameter space can be covered by existing and upcoming experiments. To discuss this, we detail a phenomenological discussion of the heavy axion properties next.

A robust consequence of the above mentioned class of models is the defining $G\tilde{G}$ coupling of the axion. For this purpose, we consider an effective Lagrangian,

with $\alpha_{i}=g_{i}^{2}/(4\pi)$ given in terms of SM gauge couplings, and $\alpha_{1}=5/3\alpha_{Y}$, in terms of the hypercharge gauge coupling. To illustrate the significance of a non-zero $c_{3}$, we will focus on two scenarios which are complementary to the case of photon or electroweak-dominance, $c_{2},c_{1}\gg c_{3}$, frequently assumed in the literature due to its testability. In more detail, we will focus on the cases of,

• •

  Gluon dominance: $c_{3}=1,c_{1},c_{2}=0$;

• •

  Codominance, $c_{1}=c_{2}=c_{3}$.

Both of the above cases are motivated from the generic Axion considerations, as well as UV considerations. In fact, these two cases match well respectively to the KSVZ Kim:1979if ; Shifman:1979if and DFSZ Dine:1981rt ; Zhitnitsky:1980tq scenario of the minimal axion theory.

These choices have an important effect on the phenomenology of such heavy axions, since for $m_{a}\gtrsim 1$ GeV, the axions predominantly decay into hadronic final states, as opposed to diphoton final states on which a significant number of searches rely.³³3For a recent study in the $aF\tilde{F}$ dominance at DUNE ND, see Ref. Brdar:2020dpr . Consequently, interesting parts of the axion parameter space open up, as we will see below. Simultaneously, a non-negligible $c_{3}$ gives rise to important axion production channels at the LHC and various proton beam dump experiments.

Below the scale of electroweak symmetry breaking, the EFT in Eq. (8) gives rise to an axion photon coupling,

with Georgi:1986df ; Bauer:2017ris ; Aloni:2018vki

To obtain $c_{\gamma}$ for $m_{a}\lesssim\Lambda_{\rm QCD}$, we have assumed the $\eta-\eta^{\prime}$ mixing angle $\sin\theta_{\eta\eta^{\prime}}=-1/3$ as in Aloni:2018vki while noting the significant uncertainty $\theta_{\eta\eta^{\prime}}\simeq-(10^{\circ}-20^{\circ})$ (see Sec. 15 of Ref. Zyla:2020zbs and references therein). Some results for more general mixing angles can be found in Refs. Ertas:2020xcc ; Gori:2020xvq . The factor of $1.92$ in Eq. (10) can be obtained after including higher order corrections diCortona:2015ldu on the leading order contribution due to quark masses $\frac{2}{3}\frac{4m_{d}+m_{u}}{m_{u}+m_{d}}\approx 2$. Importantly, we see that even in the absence of the a tree-level $c_{1},c_{2}$, the anomaly and axion-pseudoscalar meson mixing introduce a significant axion-photon coupling below $\Lambda_{\rm QCD}$.

The phenomenology will be largely dictated by the lifetime of the axion. For $m_{a}<3m_{\pi}$, the axion decays exclusively in diphoton final states with a width,

whereas above that threshold the hadronic decay modes open up Aloni:2018vki and quickly become dominant except regions near resonant mixing. We show the resulting lifetime of the axion in Fig. 3 which is used in the following to determine the reach of the DUNE ND.

The blue (red) line in Fig. 3 assumes the Gluon Dominance (Codominance) scenario, and both lines assume $f_{G}=1$ PeV. This lifetime is proportional to $f_{G}^{2}$. We note that the two scenarios are nearly identical for $m_{a}\gtrsim 1$ GeV but below $1$ GeV the above distinctions are quite important. For $m_{a}$ below 100 MeV, we can see that in the Codominance scenario, the axion has a larger lifetime as the corresponding $c_{\gamma}$ is smaller for the particular choices of $c_{i}=1$. Between 100 MeV and 1 GeV, two effects are important. One is the near resonance mixing with the SM mesons, which determines the dips in this lifetime plot. The other one is the cancellation between different meson mixings and the direct contributions to $c_{\gamma}$ in Eq. 10 for which we get peaks in the lifetime. We also note here, the mixing expansion in this equation at the meson pole regime should be regulated by the unitarity of the mixing matrix, which we neglected here in the equations but implemented effectively in the next section in our numerical computation.

Before moving on, we note that due to the $Z_{2}$ symmetry, the specific model described in Ref. (Hook:2019qoh, ) and above, has a massless mirror photon. While the axion can decay into a pair of mirror photons, for our phenomenlogical analysis below we will ignore this effect, motivated by the following reason. The mirror photon does not play an essential role in our set up and can be removed from the spectra if we do not copy the SM $U(1)_{Y}$ into the mirror sector. Instead, we can start with a common $U(1)_{X}$ under which both SM and mirror sectors are charged. As the mirror Higgs gets a VEV, the breaking $SU(2)_{W}^{\prime}\times U(1)_{X}\rightarrow U(1)_{Y}$ takes place without giving rise to a mirror photon. Since the flavor structure of both the SM and the mirror sector are the same, the effects of differential RG running of $\bar{\theta}$ between the two sectors are still suppressed.

Meson Mixing: To determine the axion production due to meson mixing, we simulate the Long Baseline Neutrino Facility (LBNF) beam as a 120 GeV proton beam colliding with a fixed target using Pythia8 with the “SoftQCD:all = on” option. For all of our simulations, we assume the total number of protons-on-target (POT) is $N_{\rm POT}=1.47\times 10^{22}$ over the course of 10 years.⁴⁴4This assumes ten years of operation at the nominal rate of $1.47\times 10^{21}$ POT/yr. The DUNE collaboration plans on upgrading its beam to a larger number of POT/yr during its operation, so our estimations should correspond to at most ten years of data collection. We find that approximately $2.89$ $\pi^{0}$, $0.33$ $\eta$, and $0.03$ $\eta^{\prime}$ are produced per POT at this beam energy.

As alluded to in Eq. (10), axions mix with the SM pseudoscalar mesons through the $G\tilde{G}$ coupling. Here we summarize the mixing angles Bauer:2017ris ; Aloni:2018vki ; Ertas:2020xcc ,

where $f_{\pi}\approx 93~\rm MeV$. In these equations, the ellipses contain $\pi-\eta$ and $\pi-\eta^{\prime}$ mixing terms, which subdominantly contribute to ALP production considered below.

Using Eqs. (12), (13) and (14), the number of axions produced from ALP-meson mixing is obtained as ⁵⁵5We ignore possible interference effects between different meson-mixing modes in this approximation.,

where

The above function $f(m_{\mathrm{meson}},m_{a})$ models the QCD production rate of mesons which decreases as one increases the meson mass, rooting from both the running strong coupling as well as the parton evolution. The power of $-1.6$ comes from fitting the $\pi^{0}$, $\eta$, and $\eta^{\prime}$ meson production rate as a function of their masses. Furthermore, we conservatively bound the function value by unity, by neglecting the possible enhancement of rate beyond the mixing calculation in the regime of $m_{a}\leq m_{\mathrm{meson}}$. Note that to our knowledge we are the first to take this further step to model the kinematic effect of masses in this regime of the axion production rate. To show the difference, we show the flux in Fig. 4 with and without this mass effect taken into account in thin and thick curves.

The axion flux with axion mass below 1 GeV at the DUNE ND then depends on Eq. (15), the ND cross-sectional area, and the acceptance fraction. In detail, for each simulated $\pi^{0}$, $\eta$, and $\eta^{\prime}$, we convert them into an axion with a weight according to the mixing angle. This conversion process keeps the energy of the SM meson the same, rescaling the magnitude of the three-momentum while maintaining their direction. The acceptance fraction is defined as the fraction of produced axions that are traveling in the direction of the DUNE ND upon production, folding in the production angular dependence in the beam dump environment. We discuss the details of the DUNE ND in our simulations in Section III.2, and we note here that the acceptance fraction for the different meson-mixing production mechanisms is $\mathcal{O}(10^{-2})$.

Fig. 4 displays the expected flux at the DUNE ND for axion production from meson mixing as a function of the mass $m_{a}$. We show the separate contributions from $\pi^{0}$, $\eta$, and $\eta^{\prime}$ as different colors, which allows us to see the different axion-meson mixing dominating for different regions of $m_{a}$. We dash the contributions for $m_{a}\gtrsim 1$ GeV, where we expect gluon-gluon fusion to serve as a better description of axion production in this environment. This flux is shown for $f_{G}=1$ PeV and scales with $f_{G}^{-2}$.

Gluon-Gluon Fusion: Above $\mathcal{O}$(GeV) masses, the direct production mode from gluon-gluon fusion could potentially dominate the contribution to the heavy axion flux at beam dump facilities. There, given the momentum exchange is above the GeV scale, the parton distribution function description is valid. The operator $\propto aG\tilde{G}$ determines the production rate.

We evaluate the production cross section convoluted with the leading order parton distribution function NNPDF Ball:2014uwa ; Hartland:2012ia with our calculation available at this link, following:

where $\mu_{F}$ and $\mu_{R}$ are the factorization and renormalization scale, and $s$ is the center of mass energy, approximately $(15~\,{\rm GeV})^{2}$. The factorization and renormalization scale is set at the heavy axion mass $\mu_{R}^{2}=\mu_{F}^{2}=m_{a}^{2}$. In Fig. 5(left), we show the inclusive production rate for a decay constant $f_{G}=1~$PeV. We used the 1-loop renormalization-group-equation running of the strong coupling constant $\alpha_{S}(Q^{2})$ embedded in NNPDF. In Fig. 5(right), we show the same rate in a linear scale, effectively zooming in around $m_{a}\approx 1$ GeV. We also show how the cross section varies when we consider factorization scales between $m_{a}^{2}$, $1/2m_{a}^{2}$ (dashed lines) and $2m_{a}^{2}$ (dotted). In both panels, we see the scale uncertainty is sizable, especially below the GeV scale. In fact, this regime is where PDFs have large uncertainties and scheme dependence. We observe with our current NNPDF choice the cross section only starts to dominate when $m_{a}\gtrsim$ 0.9 GeV, and so our results are not subject to the large uncertainty in the low mass regime. Another interesting phenomena is the cross-over of different scale choices, this is a result of the PDF evolution and the $\alpha_{S}^{2}(Q^{2})$ running from the production cross section.

In Fig. 6, we show the normalized energy distribution from the gluon-gluon fusion production for a benchmark heavy axion with a mass of 1.2 GeV. The red curves show the distribution at production, which is universal for all axion decay constants. For $m_{a}\ll\sqrt{s}$ , where $\sqrt{s}\simeq 15$ GeV for the LBNF-DUNE beam under consideration here, the distribution from gluon-gluon fusion are similar. In the same plot, we also show the differential distribution in orange, yellow and cyan for the axion accepted by the DUNE ND (accounting for decays of axions en route) for axion decay constant $f_{G}$ of 50, 100, and 200 PeV, respectively. The shift in distribution is mainly driven by the necessary boost for a heavy axion to arrive at the DUNE ND before decaying. More boost is needed for a smaller decay constant, and hence the distribution shifts towards higher energies. We discuss the detection considerations and details in Section IV.

A few other axion production modes could be significant. These include bremsstrahlung effects from the proton-proton collinear emission with resummation. The result will depend on how one treats the finite mass effect from the axion and its derivative coupling. Similarly, there can be collinear emissions from the quarks and gluons in the collision, involving model-dependent axion-quark couplings. Through the axion coupling to quarks, possible flavor-changing decays from mesons will also contribute to axion production MartinCamalich:2020dfe . Last but not least, the proton-proton collision at beam dumps will create secondary collisions from the remnants of the first collision, enlarging the number of mesons produced and hence enriching the flux for heavy axions below $\mathcal{O}\mathrm{(GeV)}$. All these effects could help improve the axion flux and, therefore, the DUNE ND sensitivities to heavy axions. We leave a detailed analysis of these different contributions with more model dependence for future studies.

We are interested in signatures of axion decay in the DUNE ND Complex. Specifically, we consider such signatures inside the liquid argon time-projection-chamber detector ArgonCube and the gaseous argon time-projection-chamber Multi-Purpose Detector (MPD). ArgonCube is situated at a distance of 574 m from the DUNE proton target and has a total active volume of $7$ m wide, $3$ m high, and $5$ m long. Fiducialization reduces this to a fiducial mass of roughly $67.2$ t Abi:2020evt ; Abi:2020kei . The MPD is situated directly downstream of ArgonCube (designed to be a spectrometer of muons and other particles that do not stop in ArgonCube) with a cylindrical volume that is roughly 5 m in diameter and 5 m in height. This corresponds to an active mass of $1$ ton. The MPD is situated inside an electromagnetic calorimeter and a magnetic field, allowing for precision measurement (and charge and particle identification) of the particles traveling through its fiducial volume. Fig. 1 provides a schematic drawing of the DUNE target and Near Detector Complex (note that many elements are removed from this figure for simplicity, including the magnetic focusing horns and a significant amount of earth between the decay volume near the target and the detector hall).

Given the dimensions of the detectors and the DUNE target/detector distance, we find that $\mathcal{O}(10^{-2})$ of the axions produced via meson mixing will travel in the direction of the near detector complex. We include axion decays inside both the liquid and gaseous detectors in our simulations, corresponding to a total decay length of roughly $10$ m. In Section IV, we discuss the various experimental signatures of this heavy axion decay and how we can reduce associated backgrounds in the two detectors.

We first focus on the case of gluon dominance discussed in Sec. II.3 where $c_{3}=1,c_{1},c_{2}=0$.⁸⁸8There would still be a non-negligible, $c_{3}$-induced photon coupling even with $c_{1}=c_{2}=0$. Combining both the meson mixing and gluon-gluon fusion production modes, we determine the parameter space for which we would expect three or more signal events in ten years of data collection at DUNE ND, the red shaded region in Fig. 10. The blue shaded regions correspond to the same ones shown in Fig. 2 based on theoretical considerations regarding the axion Quality Problem. The brown shaded regions are the same as in Fig. 2 as well. The horizontal shaded region labeled “Colored Particles” is disfavored since from a UV perspective, $aG\tilde{G}$ coupling generically originate after integrating out colored fermions with masses $\sim yf_{a}$. Requiring a maximal Yukawa coupling $y\sim 4\pi$ along with the LHC constraints on colored states to have masses above $2$ TeV Aaboud:2017nmi ; Sirunyan:2018xlo ; Sirunyan:2018rlj ; Aad:2019hjw gives the above bound.

We also include a number of existing experimental/observational constraints on this parameter space in grey. ⁹⁹9Here we update some of the astrophysical and cosmological bounds used in Ref. Hook:2019qoh by using the more recent results from Depta:2020wmr ; Ertas:2020xcc . For small mass $m_{a}\lesssim 100$ MeV, astrophysical and cosmological constraints are relevant — the region labeled “SN1987A” indicates the region of parameter space for which such axions would cool the supernova and carry away too much energy Chang:2018rso ; Ertas:2020xcc ,¹⁰¹⁰10To illustrate the uncertainty of the SN1987A bound, we show two contours corresponding to fiducial profiles used in each of Chang:2018rso and Ertas:2020xcc .¹¹¹¹11For a recent update on SN1987A bound on ALP-photon coupling, see Lucente:2020whw . whereas the region labeled “BBN + $N_{\rm eff}$” indicates where such axions would affect light-element abundances and contribute to the number of effective degrees of radiation Depta:2020wmr ; Millea:2015qra .¹²¹²12For a recent discussion on BBN constraint on ALP-lepton couplings see Ghosh:2020vti . Both the “SN1987A” bound and “BBN + $N_{\rm eff}$” are shown via dashed lines because of their associated uncertainties, see e.g. Bar:2019ifz and Depta:2020wmr , respectively. For $m_{a}$ between 1 MeV and 1 GeV, a number of searches have been performed in the contexts of both electron and proton beam dumps and corresponding bounds were discussed in Dobrich:2015jyk ; Dolan:2017osp ; Ariga:2018uku for these types of axions, as well as searches for rare meson decays Izaguirre:2016dfi ; Aloni:2018vki ; Gavela:2019wzg ; Ertas:2020xcc ; Gori:2020xvq . Furthermore, there are constraints from LHC dijet searches for $m_{a}>50$ GeV as obtained in Mariotti:2017vtv .

Other planned experiments with similar timescales as DUNE are capable of performing searches in this region of parameter space. We include some projections of these in Fig. 10 as well.¹³¹³13For comparison against other potential future experimental searches in this parameter space, see Fig. 12 in Appendix A. The FASER Ariga:2018uku (dot-dashed cyan) experiment at the LHC will be able to probe a similar mass regime as DUNE but with smaller $f_{G}$ due to its close proximity/high-energy production source. A proposed displaced decay search using the high-luminosity LHC track trigger Hook:2019qoh can probe the region encompassed by the dot-dashed brown line at heavier masses than DUNE. Finally, for $m_{a}<m_{\pi}$ , NA62 has powerful sensitivity (dot-dashed purple) via the search for $K^{+}\rightarrow\pi^{+}+a$ where $a$ is undetected Ertas:2020xcc . Comparing our DUNE projections against these other future proposals, the complementarity of these different search strategies is obvious – the combination of all of these will allow for considerable reach in the theoretically- motivated parameter space in a way that no individual experiment can accomplish on its own. DUNE will specifically be most powerful in this long-lived region of small $g_{agg}$ or large $f_{G}$, especially for $20~\text{MeV}\lesssim m_{a}\lesssim 2~\text{GeV}$.

Here we focus on the scenario where $c_{1}=c_{2}=c_{3}$ and derive the coverage for 3 events for 10 year data taking at DUNE ND, shown in red in Fig. 11. While this coverage is similar to the one in Fig. 10, around $m_{a}\sim 400$  MeV, the effective photon coupling $c_{\gamma}$, in Eq. (10) becomes small for our choices of $c_{i}$. As a result, the coverage region shifts upwards exhibiting a “peak”-like feature.

The theoretical constraints from the axion Quality Problem remain the same. The set of experimental/observational constraints are shown in grey. Since the “SN1987A” constraints are dominated by the $G\tilde{G}$ coupling, it remains the same. However, the “BBN+$N_{\rm eff}$” constraint is dominated by the $F\tilde{F}$ coupling, and hence it gets modified based on Eq. (10) as a result of non-negligible $c_{1}$ and $c_{2}$. The constraints from partially invisible kaon decays are also modified due to the non-negligible $aW\tilde{W}$ coupling Izaguirre:2016dfi ; Gavela:2019wzg ; Ertas:2020xcc . We recast the electron and proton beam dump results from Dolan:2017osp as appropriate for the present case of Codominance. Since the diphoton decay modes are now non-negligible, for higher masses $m_{a}>20$ GeV, both the diphoton and dijet searches at the LHC give relevant constraints Mariotti:2017vtv .

We also include future projections from diphoton searches at LHCb CidVidal:2018blh and HL-LHC Mariotti:2017vtv (dot-dashed green), kaon decay searches at NA62 Gavela:2019wzg ; Ertas:2020xcc (dot-dashed purple) and LHC Track Trigger proposal Hook:2019qoh (dot-dashed brown) can cover complementary regions of parameter space as before.

Some comments regarding a few omissions in Fig. 11 are in order. We expect some part of the parameter space for $100~\text{MeV}<m_{a}<1~\text{GeV}$ would be covered by the existing CHARM data Bergsma:1985qz which we have not derived for the Codominance scenario. Also, we have not derived the constraints and projections from KOTO and NA62/48 from visible kaon decays for this scenario, which would cover some parameter space for $150~\text{MeV}<m_{a}<350~\text{MeV}$ and roughly $1/f_{G}>10^{-4}~\text{GeV}^{-1}$, mostly complementary to our DUNE ND coverage. In Ref. Gori:2020xvq , such constraints were derived for the cases $G\tilde{G}$ and $W\tilde{W}$ -dominance separately. Some complementary coverage, again for roughly $1/f_{G}>10^{-4}~\text{GeV}^{-1}$ and $m_{a}\lesssim\text{few GeV}$, would also come from $B\rightarrow Ka$ processes at Belle-II similar to what is discussed in Gavela:2019wzg for the case of $W\tilde{W}$-dominance.

To summarize, similar to the case of “Gluon Dominance” above, we see that DUNE ND would provide a powerful coverage, complementary to other existing and projected constraints, especially for large $f_{G}$ and $30~\text{MeV}\lesssim m_{a}\lesssim 1~\text{GeV}$.