A Cosmic Window on the Dark Axion Portal

Axions and dark photons are common in many extensions of the Standard Model. The dark axion portal -- an axion coupling to the dark photon and photon -- can significantly modify their phenomenology. We study the cosmological constraints on the dark axion portal from Cosmic Microwave Background (CMB) bounds on the energy density of dark radiation, $\Delta N_\text{eff}$. By computing the axion-photon-dark photon collision terms and solving the Boltzmann equations including their effects, we find that light axions are generally more constrained by $\Delta N_\text{eff}$ than from supernova cooling or collider experiments. However, with dark photons at the MeV scale, a window of parameter space is opened up above the supernova limits and below the experimental exclusion, allowing for axion decay constants as low as $f_a \sim 10^4$ GeV. This region also modifies indirectly the neutrino energy density, thus relaxing the cosmological upper bound on the sum of neutrino masses. Future CMB measurements could detect a signal or close this open window on the dark axion portal.


Introduction
Standard Model (SM) extensions featuring dark sectors arise in a number of theories beyond the Standard Model (BSM).The evidence for dark matter and neutrino masses also implies the existence of a dark sector.Such a dark sector may consist of a broad spectrum of particles, in the same way that our visible sector is described by a highly non-minimal SM.This motivates studying the potential phenomenological effects of an extended dark sector for a variety of nonminimal scenarios, to anticipate ways in which they may manifest themselves in observations and experiments.
Astrophysical and cosmological observations place some of the strongest constraints on light dark sectors up to masses around the MeV scale and at weaker couplings than accessible to direct searches in experiment.Supernovae emissions in particular are sensitive to additional weakly coupled new degrees of freedom that provide an extra source of energy loss beyond neutrinos [1].Since neutrinos have been detected from the explosion of supernova SN1987A [2][3][4], supernova cooling has been used to constrain a plethora of BSM scenarios.However, there is an upper bound on supernovae exclusion limits when the couplings become strong enough to no longer enable efficient cooling, and excluding new physics in this process is subject to uncertainties in supernova modelling.Cosmological bounds on dark radiation from the Cosmic Microwave Background (CMB) are more robust.They have previously been applied to dark photons and axions, for example, in Refs.[5,6].We investigate here their implications for a nonminimal scenario consisting of two light dark degrees of freedom.In the dynamics of cosmological evolution, described by the Boltzmann equation, there may be a non-trivial interplay between the energy density in visible radiation, dark radiation, and neutrinos that lead to unconstrained regions of parameter space which would naively have been excluded.We point out such a possibility here in a BSM extension involving both an axion and a dark photon.
Axions and dark photons are two of the most widely studied dark sector candidates (see for example Refs.[7,8] for a review).Either of them could constitute dark matter, but in any case they could still be part of the dark sector.The existence of both in a non-minimal dark sector implies an axion-photon-dark photon interaction known as the dark axion portal; in an Effective Field Theory (EFT) description, all terms allowed by symmetries are a priori present.The Lagrangian term for the dark axion portal interaction is given by where a is an axion-like particle with mass m a , F µν is the field strength of the massless photon A µ , F ′ µν is the dual field strength of the dark photon A ′ µ of mass m γ ′ , and f a is the axion decay constant.The relative strengths of the various EFT interactions depend on the ultra-violet (UV) completion.Here we take the dark axion portal to be the dominant interaction between visible and dark sectors.For example, one may assume a charge conjugate symmetry in the dark sector under which the axion is odd to justify the axion-photon-dark photon coupling as the leading interaction while preventing the kinetic mixing of the dark photon to the SM.The dark axion portal of Eq. ( 1) has previously been explored in various contexts [9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26].It may also play a role in a mechanism for cosmological relaxation of the Higgs mass [20].
In this work, we compute the axion-photon-dark photon collision terms when solving the Boltzmann equations to determine ∆N eff , the additional relativistic degrees of freedom in light species other than neutrinos.The effect of the dark axion portal on the allowed dark radiation contribution from a light axion is drastically altered by the presence of a dark photon when the dark photon's mass is around the MeV scale.We find that in the excluded region f a ≲ 10 7 GeV, there is an open window of parameter space when m γ ′ is between 1 and 10 MeV where f a can be as large as 10 4 GeV.This is due to the dark photon decay after the time of Big Bang Nucleosynthesis (BBN) but before the CMB recombination era, which proceeds through the dark axion portal interaction to a massless photon and a light axion.The contribution to the visible radiation is sufficient to indirectly alter the neutrino energy density away from its value in the SM.The axion dark radiation contribution to ∆N eff can then be larger while remaining in agreement with CMB constraints from Planck data.This island of parameter space is only partly constrained by supernova cooling bounds.For stronger couplings, or equivalently smaller axion decay constants f a ≲ 10 5 GeV, supernova trapping prevents extending the supernova exclusion reach.Constraints from collider experiments, in particular CHARM [27], only enter below f a ≲ 10 4 GeV.The experimental sensitivity can be improved by the proposed SHiP experiment [28,29], though not enough to cover the mass-coupling range necessary to close the island region.There is therefore an open window of parameter space for the dark axion portal that is currently unconstrained but could be comprehensively probed by future CMB measurements such as CMB-S4 [30].
The alteration of the neutrino energy density in this window also relaxes the cosmological upper bound on the sum of neutrino masses.In particular, for certain cases this may reopen the possibility of inverted ordering which is currently excluded in the standard scenario by Planck constraints [31].A complementary probe is provided by the next generation of neutrino experiments.They are projected to reach the sensitivity necessary to cover the entire inverted ordering parameter space.Should inverted ordering by preferred, the dark axion portal has the potential to alleviate a small tension between the cosmological bound and terrestrial determinations of neutrino masses.
This paper is organised as follows: in Section 2 we review the Boltzmann equations that are used to describe the cosmological evolution of the energy densities of light relativistic species.In Section 3 we describe the results of our numerical simulation to determine ∆N eff for the case of a massive dark photon and light axion, and vice versa.Finally, in Section 4, we show how the modification of the neutrino energy density in the open window of parameter space can relax the cosmological bound on the sum of neutrino masses.In Section 5, we discuss the cosmological implications for Big Bang Nucleosynthesis (BBN).We conclude in Section 6. Expressions for the collision terms are collected in the Appendix.

Boltzmann equations
The cosmological evolution of the energy densities of the degrees of freedom present in the thermal bath of the early universe is well described by the Boltzmann equations.The energy density of relatistivic species is parametrised in terms of N eff , the effective number of neutrino species.We solve the Boltzmann equations numerically starting with a temperature around the muon mass when it is a good assumption that the universe was in a plasma state consisting of photons, electrons, positrons, and neutrinos in the SM.We include extra relativistic BSM species that can alter the cosmological evolution of the energy density and lead to an interesting cosmological signature The Boltzmann equation takes the form, where ρ i is the energy density and C i a is the collision term of a given species i for the process a.The Boltzmann equations for the SM are taken from Refs.[32][33][34][35].The additional Boltzmann equations for the dark photon and axion and the modification of the SM part in the dark axion portal model are given by where the explicit forms of collision terms are given in Appendix A. The collision terms involve five temperatures, which are different a priori, namely, T νe , T νµ = T ντ for neutrinos, T eγ = T γ = T e for the electron-photon (see e.g.Ref. [32] for related discussion), T a for the axion, and T γ ′ for the dark photon.Since we start with an initial temperature around the muon mass, only electrons were included in the Boltzmann equations and the electron neutrinos interact differently from muon and tau neutrinos.The 3-point collision term gets a contribution from the process depicted in Fig. 1.The diagrams for the leading 4-point collision terms in the coupling strength proportional to 1/f a are shown in Fig. 2.
In numerically solving the Boltzmann equations, we do not include the chemical potential (which is negligible for the electron and neutrino and vanishes for the axion, photon and dark photon) and make some assumptions regarding the statistics of the particles in the collision terms.For the process of 1 ↔ 2 + 3 with the particle 1 being the heaviest, we take into account the quantum statistics only for particle 1 and assume the Maxwell-Boltzmann distributions for particles 2 and 3, such that the integrand of the collision term is proportional to , where f i is the distribution function.For the process of 1 + 2 ↔ 3 + 4, we assume the Maxwell-Boltzmann distributions for all particles which is equivalent to setting f 1 f 2 − f 3 f 4 in the integrand of the collision terms.The above assumptions are valid at least for the heavy dark photon (or heavy axion) scenario in the mass range above 1 MeV.The 4-point collision terms can be dominant only at a high temperature where all the particles are well approximated as relativistic following the Maxwell-Boltzmann distributions.At a lower temperature, the 3-point collision term becomes dominant and the lighter particle and photon can be approximated with the Maxwell-Boltzmann distribution, whereas the quantum statistics of the heavier particle is still taken into account.New degrees of freedom in the thermal bath during Big Bang Nucleosynthesis (BBN) can modify the observed primordial abundances of helium and deuterium, for instance, through the Boltzmann equations of protons and neutrons [36]: where Note that the Hubble parameter, possibly modified due to new degrees of freedom, does not appear in Eq. 4. Since the overall Boltzmann suppression factors on both sides of Eq. 4 cancel out due to the BBN temperature being much smaller than the nucleon mass scale, they prevent terms in Eq. 4 from influencing those in Eq. 2 from which the neutrino and photon temperatures are dominantly determined.However, the modified temperatures of neutrinos and photons feed into the Boltzmann equations in Eq. 4 and may alter the cosmological evolution of primordial abundances.This will be discussed in detail in Section 5.While the collision term Γ(γ ′ + p ↔ a + p), mediated by a t-channel photon, is allowed, it will not alter the evolution of the proton abundance.The collision term will also not affect the evolution in Eq. 2 around the BBN epoch, due to the Boltzmann suppression.
A naive dimensional analysis estimate for the BBN bound was placed in Ref. [19].They find f −1 a ≲ 6.5 × 10 −7 GeV −1 for massless axions and m γ ′ < 1 MeV.The constraint is based on requiring Γ ≤ H, where Γ is the thermally averaged interaction rate around T = T BBN for e + e − → aγ ′ .This requirement effectively puts a bound on ∆N eff as it prevents axions and dark photons from being thermalised such that axions and dark photons do not contribute to the number of effective neutrinos.However, as Γ ∝ T 3 and H ∝ T 2 , the BSM species can be thermalised at higher temperatures even if they do not at BBN.In this case the cosmological constraints from N eff at the time of CMB are more relevant than BBN.In the next Section we will focus on computing the energy density of relativistic species at CMB temperatures.

Cosmological constraint from N eff
Using the Boltzmann equations described above to solve numerically for the energy density of radiation, we may now place cosmological constraints on any extra contributions from dark sectors.The effective number of neutrino species is defined as We adopt the constraint N eff = 2.99 +0.34 −0.33 at 95% CL as determined by Planck 2018 combined with polarisation, lensing, and Baryon Acoustic Oscillation (BAO) data [37].The value in the SM was recently updated to be 3.043 [38][39][40][41][42][43][44][45].
We take the axion mass to be lighter than about µeV to avoid exceeding the current dark matter relic abundance, while the dark photon mass is fixed at around the MeV scale.The opposite scenario with a light dark photon and MeV-scale axion will be considered at the end of this Section.We will not pursue the scenario where both the axion and dark photon are stable and massive enough to behave like non-relativistic particles, which would typically require an additional decay channel to avoid having too large a dark matter relic abundance.
The MeV-scale dark photon that is produced in the thermal bath when the temperature was sufficiently high will decay away at some point between BBN and CMB times.The decay proceeds through the dark axion portal to photons and axions, and thus contributes to both visible and dark radiation.In addition to the direct dark radiation contribution to N eff from the axion, the modified photon energy density means that there is also an indirect modification to the subsequent evolution of the neutrino energy density.The net effect of this joint change is to allow for a wider range of axion and neutrino energy densities at CMB recombination than would otherwise be compatible with the observed N eff .
This effect can be seen qualitatively by expressing N eff in terms of the neutrino, photon, and axion temperatures, assuming the axions to be relativistic when they decoupled, as where the first term takes into account a different neutrino-to-photon temperature ratio from the typical value (4/11) 1/3 in the SM depending on the dark photon mass and the decoupling epoch of the light axion, and the second term represents the axion dark radiation contribution with g a and g γ parametrising the number of degrees of freedom of the axion and photon respectively.The neutrino and axion temperatures consistent with observation is illustrated in the left panel of Fig. 3.We see that a modified neutrino temperature allows higher values of the axion temperature that would otherwise be disfavoured for a SM-like neutrino temperature.The points on the plot correspond to the points shown on the right panel of Fig. 3 in the plane of m γ ′ vs 1/f a .In that figure, the grey region within and below the blue solid line contours  denote the allowed N eff < 3.33 parameter space as computed from our numerical simulation.For an axion in thermal equilibrium, the axion temperature T a will be the same as the photon radiation temperature T γ .A smaller axion temperature can only be achieved by decoupling the axion from the thermal bath earlier, which happens for sufficiently small values of 1/f a .The dark photon in that case would also not be produced appreciably (for instance, see the bottom panels of Fig. 4 that will be discussed later).This leads to the horizontal blue exclusion line below which the parameter space is unconstrained.On the other hand, for the situation with a smaller T ν /T γ than the typical value in the SM, a larger axion-to-photon temperature ratio can be more easily achievable.This happens automatically when the dark photon with mass near the BBN scale decays to the axion and photon around or after neutrino decoupling.The entropy of the dark photon is dumped mostly into the photon sector, thus decreasing the neutrino-to-photon temperature ratio.An increased T a /T γ is then compensated by a smaller T ν /T γ in the overall contribution to N eff .The decoupling of the axion from the thermal bath can be delayed, allowing for stronger couplings with a smaller value of f a .In this situation, we expect to see a nontrivial exclusion curve near the BBN scale in the plane (m γ ′ , 1/f a ), as is confirmed by the island region in the right panel of Fig. 3.The horizontal exclusion strip around roughly 1/f a ∼ 10 −7 GeV −1 that separates the two allowed regions occurs due to an intermediate regime where the dark photon is no longer sufficiently abundant to allow for a large enough modification of the neutrino temperature but still contributes an appreciable amount of axion dark radiation.For weaker couplings, at larger GeV.The green band corresponds to the allowed region, N eff = 2.99 +0. 34 −0.33 from Planck 2018 [37].Right: the evolution of the i-species-to-photon temperature ratios for the same benchmark scenarios.The grey horizontal line corresponds to T ν /T γ = (4/11) 1/3 .f a , neither the dark photon nor axion are produced significantly enough to affect CMB bounds.We note that the size of the horizontal exclusion strip depends mildly on the assumption of the reheating temperature which determines the freeze-in abundance of dark photons in the thermal bath.In our numerical simulation of the Boltzmann equations, the initial temperature is set to be around a few hundreds of MeV as the evolution for heavier SM particles than the electron are not included.We checked that increasing the initial temperature to a higher value beyond the valid region for our numerical simulation widens slightly the horizontal exclusion band around 1/f a ∼ 10 −7 GeV −1 as more dark photons freeze-in to the thermal bath, but this widening tapers off.If we, instead, decrease the initial temperature down to roughly 50 MeV, the horizontal exclusion band shifts upwards slightly.As illustrated in Fig. 5, the horizontal exclusion band moves up with decreasing initial temperature while the upper part of the unconstrained island of parameter space remains intact.This horizontal band therefore has a mild dependence on the initial reheating temperature.We discuss further the underlying physics with supplementary plots in Appendix C and Fig. 13.
We can see more quantitatively the aspects discussed above by following the cosmological evolution of the individual energy densities and temperature ratios.Three benchmark points with f a = 10 5 , 10 6 , and 10 7 GeV and m γ ′ = 4 MeV, corresponding to the black triangle, circle and diamond points of Fig. 3, are illustrated in Fig. 4 (similar plots for varying m γ ′ with fixed f a are illustrated in Fig. 11 in Appendix C).The left and right panels show the evolution of supernova cooling (blue) [19], CHARM (green) [27], and BaBar (orange) [46].The projected sensitivity from SHiP [28,29] is denoted by a dashed red line.
the energy densities for each species in terms of their N eff and their temperature ratio with respect to the photon temperature, respectively.In the left panels, the horizontal green band represents the range of N eff that is compatible with observation, and the solid black line is the total value of N eff that we computed.The dashed black and dot-dashed black lines represent the N eff contributions from neutrinos and axions, respectively, and the solid red line shows the evolution of the dark photon energy density.We see that the dark photon decay affects indirectly the neutrino energy density that evolves to lie outside the green band, but the dark radiation contribution compensates to allow the total N eff to remain compatible with CMB observations.In the right panels, the horizontal grey lines correspond to the expected value of the neutrino to photon temperature ratio in the SM, T ν /T γ = (4/11) 1/3 .The solid black line is the electron neutrino temperature, while the dashed and dot-dashed lines correspond to the other neutrino species' temperature and the axion temperature, respectively, all relative to the photon temperature.The enhancement in the dark radiation from the axion and the deviation of the neutrino temperature from the SM case is clearly visible.The top and middle panels confirm our picture of the physics leading to the island region.The bottom panels illustrate the situation when the dark photon freezes in at weaker couplings, which leads to the horizontal plateau in the right plot of Fig. 3.In Fig. 6, our allowed (white shaded) and excluded (grey shaded) regions in the plane (m γ ′ , 1/f a ) is overlaid with other constraints from supernova cooling and collider experiments as presented in Refs.[19,24,25,47].The constraint on 1/f a ∼ 10 −7 GeV −1 corresponding to the horizontal plateau is stronger than the conservative BBN bound derived in Ref. [19], and reaches the lower end of the blue-shaded constraints from supernova emission of Ref. [19].The BaBar [46] and CHARM experiments [27] are shaded in orange and green respectively.We see that they are not able to probe the allowed parameter space in the island region that lies above the supernova constraints.The proposed SHiP experiment [28,29], on the other hand, could have the necessary sensitivity to cover some of the parameter space, as shown by the dashed red line.
Finally, we show in Fig. 7 the plots corresponding to those in Fig. 3 but for the opposite scenario of a heavy MeV-scale axion with a light dark photon, with the same axes as in Fig. 3. Since the collision terms involved are identical, the only difference is in the number of degrees of freedom.Replacing g a with g γ ′ and T a with T γ ′ in the expression for N eff in Eq. 6, we see that there is a factor of two enhancement from the relativistic dark photon having two degrees of freedom instead of one as for the axion.The modified neutrino temperature is then no longer sufficient to compensate for the dark radiation contribution, as shown by the triangle and square black points on the left plot that are in the excluded white region in the right plot.In this case there is no island feature in the exclusion plane.

Relaxing the cosmological bound on neutrino masses
In the island region identified in the previous Section 3, the modified neutrino energy density leads to a relaxed CMB bound on neutrino masses compared to the standard scenario with no extra dark radiation.This can be seen by the expression for the neutrino energy density, ρ ν = m ν n ν , when neutrinos become non-relativistic, where the sum of neutrino masses can be factorised out if the number density is universal (which is a good assumption since the difference between the muon and tau neutrino energy densities and the electron neutrino energy density is negligible according to Fig. 4).Using CMB measurements, the neutrino abundance with respect to the critical density ρ 0 crit is approximately given by where h is the reduced Hubble parameter.The Planck 2018 data together with BAO, polarisation, and lensing leads to a strong upper bound on the sum of neutrino masses [37], m ν < 0.12 eV at 95% CL.A more recent combination from Ref. [31] leads to a stronger upper bound, m ν < 0.09 eV at 95% CL.Taking the latter constraint, a reduced neutrino temperature relative to the photon temperature, T ν /T γ , will then roughly relax the bound on neutrino masses by a factor of n SM ν /n ν , where n ν ∝ (T ν /T γ ) 3 and n SM ν is the neutrino number density in the SM.In the near future, data from the ground-based DESI (Dark Energy Spectroscopic Instrument) experiment [48] and the EUCLID satellite [49] are expected to reach a precision for the sum of neutrino masses up to 0.02 eV at the 1σ level.
Neutrino oscillation experiments [50][51][52] provide a complementary probe of neutrino masses.They constrain the difference in squared neutrino masses, ∆m 2 ν , and put terrestrial constraints on m ν as a function of the lightest neutrino mass m 0 .These are illustrated in Fig. 8, where the value of m ν is constrained to be inside the red and blue coloured band for the normal and inverted neutrino mass orderings respectively.We see that the sum of neutrino masses must be at least larger than 0.06 eV (0.10 eV) for the normal (inverted) ordering.The cosmological bounds from the CMB in the standard ΛCDM model are shown by the solid black horizontal line, which excludes the grey region above corresponding to m ν < 0.09 eV [31].The inverted ordering would appear to be ruled out by the CMB, which would lead to a tension in this scenario between terrestrial and cosmological observations.However, for a conservative benchmark, we illustrate how the CMB constraint can be relaxed in the presence of a dark axion portal with m γ ′ = 3 MeV and f −1 a = 1.6 × 10 −5 GeV −1 .In this case the neutrino to photon temperature ratio falls down to T ν /T γ ≃ 0.675.The corresponding exclusion is shown by the purple region above the horizon dashed line in Fig. 8, where the inverted ordering is now compatible with this cosmology for m 0 ≲ 0.005 eV.
The mechanism pointed out here can only mildly relax the cosmological bound from the CMB.As N eff bounds get stronger due to data from DESI and EUCLID, and in the future from proposed instruments such as CMB-S4 [30], the constraint could be further relaxed in other  is the allowed sum of the neutrino masses as a function of the lightest neutrino mass m 0 for the normal (inverted) ordering.The gray region is excluded at 95 % CL from the most recent cosmological data [31].The exclusion region is relaxed to the purple region for our benchmark point of m γ ′ = 3 MeV and f −1 a = 1.6 × 10 −5 GeV −1 in the dark axion portal model.
approaches such as those proposed in Refs.[53,54].Since they directly affect the neutrino sector, larger modifications of the neutrino temperature are possible.Other possibilities include using BSM neutrino decays [55][56][57][58][59] or other exotic modifications such as time-dependent neutrino masses or non-standard neutrino momentum distributions [60][61][62][63][64][65].Some variation of these scenarios may naturally be combined with the dark axion portal in a more UV-complete dark sector model.We leave the investigation of these possibilities to future studies.

Cosmological constraint from BBN
Consider first a simplified setting where only the neutrino decoupling temperature is modified, assuming the standard photon temperature.Then consistency with the observed primordial abundance of deuterium and helium requires a neutrino decoupling temperature T dec.ν > 0.6 and 0.3 MeV, respectively, at the 1σ level [67].This bound was derived by solving for the neutrino decoupling temperature with a modified Fermi constant G ′ F without adding any new BSM species.The modified G ′ F affects the evolution of all species including the primordial abundances.The neutrino temperature in the SM is T dec.
ν, SM ∼ 1.9 MeV, which grows with smaller G ′ F and reduces with larger G ′ F , whereas the ballpark of T dec.ν in our scenario is O(MeV) and the Fermi constant in Eq. 4 is the same as the SM.In this simplified case the primordial abundance is not expected to put any constraint on our scenario.
However, the dark photon decaying into axions and photons around the BBN temperature injects some energy into the thermal bath of the visible sector, thus modifying the photon temperature.We calculate the primordial Helium abundance by using the exact evolutions of the neutrino and photon temperatures obtained from the Boltzmann equations in Eq. 2 together with the Boltzmann equations of the protons and neutrons in Eq. 4. Following Ref. [32] (see Appendix A.4), the Helium abundance Y p was estimated as Y p ≃ 2X n | T D with T D = 0.07 MeV and we picked some representative benchmark points within the currently unconstrained parameter space.T D is the temperature at which Deuterium is no longer dissociated by photons [68,69] and thus all neutrons at that time are expected to form Helium.The resulting primordial Helium abundance Y p from our calculation is illustrated in Fig. 9 as a function of f a for fixed m γ ′ = 4 MeV (left panel), and varying m γ ′ with fixed f a = 10 6 GeV (right panel).The f a value around the bump in the left panel of Fig. 9 matches to the exclusion band in Figs. 3 and 6.We see that the benchmark points of our parameter space in Figs. 3 and 6 are compatible with the constraint from the Helium abundance, denoted by the horizontal grey band [66].The exact value of T D ≃ 0.07 MeV and the resulting primordial helium abundance will depend on a more involved numerical simulation that is beyond the scope of this work.We illustrate the sensitivity of Y p to T D in Fig. 10 for completeness.We do not expect the Deuterium abundance to lead to stronger constraints in our scenario as the evolution of the energy density in Fig. 4 and 11 indicates that the dark photon decay occurs around the temperature of 1 MeV and the photon injection continues until N eff becomes saturated around T γ ∼ 0.1 MeV.The resulting neutrino-to-photon temperature may then affect both N eff and the Deuterium abundance in a correlated way, and a large portion of the currently unconstrained parameter space may survive when the deviation of N eff from the SM value is small.Our study motivates a more comprehensive numerical analysis for the Deuterium abundance that we leave for future work.

Conclusion
The dark axion portal is an interaction between the axion, photon and dark photon.It is motivated by non-minimality of the dark sector, where axions and dark photons have long been considered potential SM extensions in many BSM scenarios.If both are present, a dark axion portal interaction term is mandated by EFT and may significantly alter the phenomenology.We computed the Boltzmann equations for the evolution of the thermal bath in the early universe including the axion and dark photon interacting dominantly with the visible sector through the dark axion portal.The dark photon is taken to be at the MeV scale, with the axion much lighter.We found that the cosmological constraints from the number of relativistic degrees of freedom, N eff , in CMB observations are generally stronger than astrophysical and collider bounds, and identified an island region of parameter space that is currently unconstrained.The weakening of the constraints in that region is due to the neutrino energy density being indirectly modified by the dark photon decays, with axionic dark radiation compensating to maintain an N eff value compatible with observation.This also has the effect of relaxing the cosmological bound on the sum of neutrino masses which no longer excludes the inverted ordering scenario that is still compatible with terrestrial neutrino mass constraints.
Future experiments such as SHiP could probe a part of this open window on the dark axion portal at relatively strong coupling, or equivalently for small decay constants f a ∼ 10 4 GeV.Better CMB observations will be necessary to cover the entire island region parameter space and further improve the sensitivity to larger decay constants.Should a potential BSM signal emerge in cosmological N eff measurements or in terrestrial neutrino mass determinations, there may be some non-trivial interplay between the two that could give us a handle on the physics of the dark sector.

A Collision terms
In this section, we derive general formulae for collision terms that we used in our analysis.While, in principle, the multi-dimensional numerical integration of collision terms should be possible, the overall performance is improved if the number of integration variables are reduced before the numerical integration.To simplify the discussion below, we distinguish two types of collision terms, denoted by C and C, contributing to the evolution of the number density and energy density, respectively.Recall that the energy density is defined as where g is a degree of freedom of the particle and f is the distribution function of a given species.
The collision term C is obtained by convoluting C with the energy E of the corresponding species.

A.1 Three-point collision term
The 3-point collision term is relevant for γ ′ ↔ aγ and the one in the evolution of the number density is given by, for the 1 ↔ 2 + 3 process, Assuming the Maxwell-Boltzmann distributions for the decayed particles 2 and 3 and the Bose-Einstein statistics for the particle 1, and the matrix element M can be factored out of integrals, the above expression simplifies to where The term in Eq. ( 11) should be convoluted with the energies E 1 , E 2 , and E 3 = E 1 − E 2 for particles 1, 2, and 3, respectively, to convert to collision terms in the Boltzmann equation for the energy density:

A.2 Four-point collision term
The 4-point collision term inducing the evolution of the number density for 1 + 2 ↔ 3 + 4 is We reduce the number of integration variables for a better analytic understanding and for better performance of the numerical evaluation.While it is a difficult task for the general situation taking into account quantum statistics for all particles, as was argued in Section 2, we assume the Maxwell-Boltzmann distributions for all particles.
where η is the polar angle of ⃗ p 2 with respect to ⃗ p 1 and θ and ϕ are polar and azimuthal angles of ⃗ p 3 with respect to ⃗ p 1 .The factor 4π (2π) accounts for the freedom in choosing the orientation of ⃗ p 1 (⃗ p 2 relative to ⃗ p 1 ).Since M 4 is given in terms of the Mandelstam variables, we convert the integrations over η and θ to those in terms of s and t whereas the integration over ϕ is removed using the delta function in Eq. (15).
where the upper and lower bounds for the integrations over t, s, and E 3 are given by, respectively, The collision term for the energy density of ith particle is given by In the evaluation of C 1 4 and C 2 4 , the integration over E 3 can be removed by noting that

B Collision terms for the massive dark photon and light axion case
We present here explicit forms of collision terms for the scenario with an MeV-scale axion and light photon, m γ ′ ≫ m a .

B.1 Dark photon decay
The collision term in the evolution of the energy density for γ ′ ↔ aγ (1 ↔ 2 + 3) is copied from the formula in Eq. ( 13): where the amplitude squared (summed over polarizations of initial and final states) and the upper/lower bound of the integration over E a are given by For each collision term, the integration over E a can be done analytically.We obtain the following collision terms for the energy density, Note that these collision terms assume the Maxwell-Boltzmann distributions for a and γ and take into account the quantum statistics for the dark photon.

B.2 s-channel process
The amplitude squared summed over spins and polarizations for the electron pair creation and annihilation process γ ′ a ↔ e + e − (1 + 2 ↔ 3 + 4) mediated by the photon is given by Assuming the Maxwell-Boltzmann distributions for all particles, we evaluate the collision term for the energy density using the formula in Eq. ( 18).The term f 3 f 4 , or f e + f e − , in Eq. ( 18) can be expressed in terms of E γ ′ and E a , The four-point collision terms for the dark photon and axion, C γ ′ γ ′ a↔e + e − and C a γ ′ a↔e + e − , respectively, have the same form except the convoluted energy in the integrand for the energy density, and the one for the electron pair is obtained using energy conservation, E γ ′ + E a = E e + + E e − .They are given by where the upper and lower bound of s are given by

B.3 t-channel process
The squared matrix element for t-channel process e ± γ ′ ↔ e ± a is given by While the simplified analytic expression of the collision term is not available, it can be approximated for the massless axion (m a = 0) in the limit of the massless electron, or m e → 0. The residual dependence on m e in logarithmic terms is to regulate the infrared divergence.They are given by C C e e ± γ ′ ↔e ± a = C γ ′ e ± γ ′ ↔e ± a − C a e ± γ ′ ↔e ± a .GeV in the heavy axion and light dark photon scenario.The green band corresponds to the allowed region, N eff = 2.99 +0. 34 −0.33 from Planck 2018 [37].Right: the evolution of the i-species-tophoton temperature ratios for the same benchmark scenarios.The gray horizontal line corresponds to T ν /T γ = (4/11) 1/3 .In the bottom panels, curves for the dark photon and axion are not visible due to their tiny contributions.GeV −1 in heavy dark photon and light axion scenario.The green band corresponds to the allowed region, N eff = 2.99 +0. 34 −0.33 from Planck 2018 [37].Right: the evolution of the i-species-tophoton temperature ratios for the same benchmark scenarios.The gray horizontal line corresponds to T ν /T γ = (4/11) 1/3 .

Figure 3 :
Figure3: Left: Allowed N eff band for the neutrino-to-photon and axion-to-photon temperature ratios assuming that the decoupling of a particle is instantaneous for the heavy dark photon and light axion scenario.Right: Exclusion plot in the plane (m γ ′ , 1/f a ) showing the contours of N eff values.The dark grey region inside the island bounded by the blue solid line and below the horizontal solid blue line with N eff < 3.33 are allowed.

1 )Figure 5 :
Figure 5: Exclusion plot in the plane (m γ ′ , 1/f a ) showing the contours of N eff values from simulations with different initial temperatures.The dark grey region inside the island bounded by the blue solid line and below the horizontal solid blue line with N eff < 3.33 are same as Fig. 3 and they were obtained with T ini = 200 MeV.Solid (dashed) red lines correspond to the allowed regions with N eff < 3.33 obtained with T ini = 100 (50) MeV.

Figure 6 :
Figure 6: Our result showing the region excluded by ∆N eff in grey, overlaid with the constraints from

Figure 7 :
Figure 7: Constrained regions for the scenario with a massive MeV-scale axion and lighter dark photon.Left: Allowed N eff band for the neutrino and dark photon temperatures relative to the photon temperature, assuming instantaneous decoupling.Right: Exclusion plot in the (m a , 1/f a ) plane showing contours of N eff values.The region below the blue-coloured line with N eff < 3.33 is allowed.

Figure 8 :
Figure 8: Cosmological and terrestrial bounds on neutrino masses.The red (blue) colored band

Figure 9 :
Figure 9:The Helium abundance Y p as a function of the interaction strength f a for m γ ′ = 4 MeV (left) and as a function of the dark photon mass m γ ′ for f a = 10 6 GeV (right) in the heavy dark photon scenario.T D = 0.07 MeV was chosen in both plots.The grey band corresponds to Y p = 0.245 ± 0.006 at 2σ[66].

= 10 6 Figure 10 :
Figure 10: The Helium abundance Y p as a function of T D for various interaction strength f a with m γ ′ = 4 MeV (left) and various dark photon mass m γ ′ with f a = 10 6 GeV (right) in the heavy dark photon scenario.The gray band corresponds to Y p = 0.245 ± 0.006 at 2σ [66].

Figure 12 :
Figure 12: Left: evolutions of N eff ≡ N ν eff + N γ ′eff and individual contributions for m a = 4 MeV and f a = 10 6 , 10 7 , 10 8 GeV in the heavy axion and light dark photon scenario.The green band corresponds to the allowed region, N eff = 2.99 +0.34  −0.33 from Planck 2018[37].Right: the evolution of the i-species-tophoton temperature ratios for the same benchmark scenarios.The gray horizontal line corresponds to T ν /T γ = (4/11) 1/3 .In the bottom panels, curves for the dark photon and axion are not visible due to their tiny contributions.