Hidden monopole dark matter via axion portal and its implications for direct detection searches, beam-dump experiments, and the H0 tension

Hidden monopole is a plausible dark matter candidate due to its stability, but its direct experimental search is extremely difficult due to feeble interactions with the standard model particles in the minimal form. Then, we introduce an axion, a, connecting the hidden monopole and the standard model particles and examine the current limits and future prospects of direct dark matter searches and beam-dump experiments. We find two parameter regions around ma = O\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{O} $$\end{document}(10) MeV, fa = O\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{O} $$\end{document}(105) GeV and ma = O\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{O} $$\end{document}(100) MeV, fa = O\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{O} $$\end{document}(104) GeV where monopole dark matter and the axion are respectively within the reach of the future experiments such as PICO-500 and SHiP. We also note that the hidden photons mainly produced by the axion decay contribute to dark radiation with ∆Neff ≃ 0.6 which can relax the H0 tension.


Introduction
Dark matter (DM) is ubiquitous in the universe, making up about a quarter of the present energy density. It is also known to play a crucial role in structure formation. However, what DM is composed of remains one of the great mysteries in particle physics and cosmology.
One of the peculiar properties of DM is its stability. Its lifetime must be at least an order of magnitude longer than the current age of the universe [1,2]. The stability of DM may be ensured by symmetry. For instance, in the WIMP scenario, its stability is usually achieved by imposing Z 2 symmetries such as R-parity and KK-parity. It is also possible that the stability of DM is due to an unbroken gauge symmetry (see ref. [3] and references therein). One such example is a hidden magnetic monopole, which is the main focus of this paper.
A monopole is a topological defect which appears associated with the spontaneous breaking of a gauge symmetry G down to H with a non-trivial π 2 (G/H). The simplest example is the 't Hooft-Polyakov monopole [4,5] with G = SU (2) and H = U(1), and it is the magnetic charge under the U(1) gauge symmetry that makes the monopole absolutely stable. If there is a hidden sector which includes such a gauge symmetry and its subsequent spontaneous breaking allows the monopole solution, then hidden monopoles are necessarily produced at the phase transition, and they can contribute to DM. 1 The hidden monopole JHEP01(2020)185 DM was studied in detail in refs. [6][7][8]. Through the Witten effect [9], the monopole DM can suppress the abundance and the isocurvature perturbations of the QCD axion [10][11][12][13].
The monopole is a promising candidate for DM, but it is formidable to search directly by experiments because its interactions to the standard model (SM) particles are feeble. In fact, in the minimal setup, it is completely decoupled from the SM sector at a renormalizable level. The setup can be extended to couple the monopole to the SM sector via the Higgs, vector, or axion portal coupling. In fact, the Higgs portal coupling was assumed in refs. [7,8] in order to keep the hidden sector in thermal equilibrium with the SM plasma. The monopole can interact with the nucleon via the Higgs portal, but the expected cross section is much smaller than any direct DM search experiments in the foreseeable future. On the other hand, the vector portal coupling induces a small electric charge for the hidden monopole. The mini-charged DM is strictly constrained by various experiments and observations [14]. We will come back to this possibility later in this paper. Our main interest is the axion portal coupling. As described below, the monopole DM scatters off nucleons with a cross section within the range of future experiments such as PICO-500 for a certain parameter region. Furthermore, axions can be probed by beam-dump experiments such as SHiP. As we shall see later, we find two parameter regions with the axion mass The hidden sector contains a massless hidden photon which contributes to dark radiation. The amount of dark radiation is customarily expressed in terms of the extra effective neutrino species, ∆N eff . Due to the axion portal coupling, the hidden photon abundance is mainly determined by the axion decay when the cosmic temperature becomes comparable to the axion mass. In particular, for the above parameter space where both monopole DM and the axion are within the experimental reach, the amount of dark radiation is predicted to be ∆N eff 0.6. We note that such an amount of dark radiation can significantly relax the H 0 tension [15].
The rest of this paper is organized as follows. In section 2, we briefly review the monopole DM and the Witten effect. In section 3, we introduce the axion portal coupling and solve for the axion profile around a monopole. Using the axion profile, we estimate the scattering cross-section of the hidden monopole with nucleons and discuss its implications for future experiments. In section 4, we discuss the current limits and future prospects of various experiments including direct DM search experiments and the beam-dump experiments. The last section is devoted to discussion and conclusions.

't Hooft-Polyakov monopoles
't Hooft and Polyakov showed that a magnetic monopole can arise when a non-Abelian gauge symmetry is spontaneously broken via the Higgs mechanism [4,5]. The simplest possibility is an SU (2) Table 1. The masses, electric and magnetic charges of the particles in the hidden sector, where α H = e 2 H /(4π) is the hidden fine structure constant. Here we take the Bogomol'nyi limit for the hidden monopole mass and its hidden electric charge is due to the Witten effect.
isovector scalar field φ = (φ 1 , φ 2 , φ 3 ). The Lagrangian density is given by H is the field strength tensor of the SU(2) H gauge fields with e H > 0 being the hidden gauge coupling constant, and D µ φ = ∂ µ φ + e H A µ H × φ is the covariant derivative. 2 The potential of the isovector scalar field is where λ φ > 0 is the quartic coupling. The potential is minimized at φ = υ H , where SU(2) H is broken down to an unbroken U(1) H . Now if we take φ = (0, 0, υ H ), the unbroken U(1) H symmetry corresponds to rotations about the 3-axis, and the hidden photon is identified with A µ H3 . Expanding the Lagrangian density around the vacuum state, we then obtain a massless hidden photon γ H , two massive hidden gauge bosons W ± H , and a massive neutral Higgs field ϕ. In addition, this theory allows a static field configuration with finite energy due to the nontrivial π 2 (SU(2) H /U(1) H ) = Z. Such a topological soliton with nonzero winding number possesses a magnetic charge, and is therefore identified with the magnetic monopole M. The hidden magnetic charge of the monopole with a unit winding number is Q M = 4π/e H .
We summarize their masses and electric and magnetic charges in table 1. In particular, the mass of the hidden monopole is bounded from below by the Bogomol'nyi bound [16], where F is a monotonically increasing function with F (0) = 1 and F (∞) 1.787 [17]. Throughout this paper, we will adopt the hidden monopole mass saturating the Bogomol'nyi limit (λ φ → 0) for simplicity following ref. [8]. The electric charge of the hidden monopole due to the Witten effect will be explained in the next subsection.

The Witten effect
To see how the hidden monopole acquires a nonzero electric charge, let us briefly review the consequence of adding the following P-and CP-violating θ H -term in the hidden sector, is the field strength tensor of the hidden photon, and F µν H is its dual tensor. We will promote θ H to a dynamical variable, the axion, in the next section. Since the θ H -term can be rewritten as a total derivative, it does not affect the equations of motion, and one usually dismisses such a term. However, in the presence of the monopole, it has a physical effect. In other words, the monopole acquires an electric charge proportional to θ H , in addition to the magnetic charge. This is known as the Witten effect [9]. As a result, the monoploe becomes a dyon [18].
To derive the explicit form of the electric charge of the hidden monopole, let us start with the equations of motion of the hidden electromagnetism, is the magnetic 4-current density. Noting that the electric and magnetic fields appear in the field strength as F j0 H = E j H and F jk H = − jkl B l H , we obtain the following modified Maxwell's equations, Suppose that there exists a monopole with the magnetic charge Q M sitting at the origin. When measured far from the origin, the magnetic field induced by the monopole is where r is the distance away from the origin andr denotes the unit vector along the radial direction. Here we assume that the core radius of the monopole is negligibly small compared to r. Using (2.5), one can see that a nonzero hidden electric field is also generated around the monopole at a sufficiently large r, Thus, the hidden monopole acquires a nonzero hidden electric charge (see table 1). Note that anti-monopoles have the opposite hidden electric charges. One can extend this result by promoting θ H to an axion field and show that the electric charge of the monopole is similarly induced by the Witten effect [19]. The difference is that the electric charge spreads out in space, and the electric charge contained inside a sphere with a radius r is given by Q E (r) = e H θ(r)/(2π). Thus, the electric charge distribution depends on the axion field configuration around the monopole which will be studied in section 3.

Relic abundance
In the 't Hooft-Polyakov monopole model, both massive gauge bosons and monopoles are stable because they carry electric and/or magnetic charges under unbroken U(1) H symmetry and there are no other lighter charged particles. Once produced, therefore, they will contribute to the DM.
The hidden sector is assumed to be in thermal equilibrium with the SM sector through some portal interactions. Then, the massive gauge bosons are in thermal equilibrium and their abundance is fixed when the annihilation processes such as W + H W − H ↔ γ H γ H , ϕϕ are decoupled as in the usual freeze-out scenario. On the other hand, monopoles are produced at the phase transition, and their abundance depends on the order of the phase transition.
Throughout this paper, we assume that the universe experiences a second-order phase transition at the critical temperature T c so that the monopole abundance is determined by the Kibble-Zurek mechanism [20,21]. 3 At high temperatures, T T c , the SU(2) H symmetry is unbroken, and the vacuum expectation value of the isovector scalar field vanishes everywhere in the universe. When T T c , the SU(2) H symmetry is spontaneously broken as the value of φ changes from 0 to υ H . Then, a characteristic domain structure appears, where the direction of φ in each domain is randomly chosen. The size of each domain is determined by the correlation length, ζ. Toward the phase transition, both the correlation length and the relaxation time scale increase proportional to |T − T c | −ν with ν being the critical exponent. The classical value of the critical exponent is ν = 1/2. The increase can be understood by noting that the dynamical time scale is determined by the inverse of the effective mass of φ, which becomes smaller and approaches zero toward the critical point. At a certain point, the relaxation time becomes equal to the time left for the universe to reach the critical point, beyond which the typical size of the domains is frozen. Then, approximately one monopole is produced in the volume of order ζ 3 . The monopole abundance from the Kibble-Zurek mechanism is thus estimated as [6,21] where h 0.7 is the reduced Hubble parameter, and M Pl 1.22 × 10 19 GeV is the Planck mass. Notice that the monopole-anti-monopole annihilation after the phase transition is taken into account in deriving (2.9).
then the hidden monopole mass should be about several hundred TeV in order to explain the observed DM density [6]. The monopole DM has sizable self-interactions because it has both magnetic and electric charges. Although there are various observational constraints on self-interacting DM, such heavy monopole DM is considered to escape the restriction [3,8]. In our numerical calculation, we will pick up α H 0.73, m M 216 TeV, and ν = 0.5 as a benchmark point, for which around 35% of the DM is contributed by the hidden monopole [8]. Note that, although the precise fraction of the monopole DM is subject to various uncertainties such JHEP01(2020)185 as the critical exponent and the presence of hidden fermions that are necessary for the interaction of the axion to be introduced in the next section, they do not change our results by more than an order of magnitude at least in the parameter region of our interest.

Spin-dependent monopole-nucleon interactions through the axion portal
In this section, we will demonstrate how to calculate the spin-dependent cross-section of the hidden monopole DM scattering off a nucleon by introducing the axion portal. Such calculation may look intractable as the monopole is a composite object and it is not clear how to express its interaction with the axion and the SM particles in the Lagrangian. In the following, we derive the axion configuration around the monopole due to the Witten effect. In the case of the massless axion, it was obtained in ref. [19]. Here we extend their result by including the axion mass, and use the obtained axion configuration to estimate the cross section. The inclusion of the axion mass is crucial because the tail of the axion profile gives the dominant contribution to the scattering cross section. If the axion has Yukawa-like interactions with the SM fermions, the monopole surrounded by the axion can scatter off a nucleon. Then, its leading-order cross section can be estimated by treating the axion configuration as an external field. This enables us to study the current limits and future prospects of various direct DM search experiments.

Axion portal
To find out the axion profile around the hidden monopole, let us consider the following Lagrangian density, where a = a(r) is the axion field with f H being the decay constant of the axion. If the shift symmetry of the axion is linearly realized as global U(1) PQ symmetry in the UV completion, it resides in the phase of a complex scalar, S. Then, the decay constant is considered to be of order the vacuum expectation value of S in a simple setup. Also, the axion coupling to the hidden photons is generated if S has a coupling to hidden fermions charged under SU(2) H . The mass of the hidden fermions is also expected to be of the same order. 4 In eq. (3.1), we have introduced the axion mass term. This may originate from another hidden gauge theory SU(N) H with gauge field strength F H , in which θ 0 is the coefficient of the theta term, θ 0 F H F H , and the mass term is obtained by expanding the potential around θ 0 . The expansion may break down when θ 0 1 as a/f H changes its value more than O(1) around the monopole. We, however, neglect the possible deviation from the quadratic potential because the calculation of the scattering cross section is essentially determined by the tail of the axion configuration where the axion potential can be well approximated 4 The massive gauge boson W ± H may decay into those fermions if kinematically allowed. For the parameter region of our interest, the mass of hidden fermions is not very different from W ± H , and so, this does not change our result significantly.

JHEP01(2020)185
by the quadratic form. Alternatively, one can simply break the shift symmetry of the axion just by introducing such a mass term. Now, let us place a monopole at the origin and derive the static field configuration of the axion around it. We suppose that, in the absence of monopoles, or at a sufficiently large distance from monopoles, the axion is stabilized at θ 0 . Due to the Witten effect, the monopole is surrounded by an electric charge, and the total amount of the electric charge is given by e H θ 0 /2π, which is independent of the axion decay constant. The electric field becomes stronger as the axion approaches the monopole. As a result, the axion field value becomes close to zero to suppress the energy of the electric field. Thus, the axion configuration is determined by the balance between the gradient energy of the axion and the energy stored in the electric field around the monopole.
To derive the precise axion distribution around the monopole, one needs to minimize the Hamiltonian of the axion-monopole system given by where we have defined a dimensionless axion field, θ ≡ a/f H , and we have used eqs. (2.7) and (3.1) and omitted the kinetic energy in the second line as we are interested in the static field configuration. We also drop the energy density of the hidden magnetic field since it does not contribute to the axion profile. 5 We have assumed spherical symmetry of the solution in the second line, and r 0 is defined as The lower end of the integration is set by the core radius of the hidden monopole, r c . In the case of 't Hooft-Polyakov monopole, it is given by r c m −1 W H , inside which the original SU(2) H gauge symmetry is restored. In the following analysis, we will assume that the core radius is negligibly small. This is because the axion profile outside the core is considered to be insensitive to r c , since the electric field outside is simply determined by the Gauss's law, and it does not depend on the core radius. Thus, in practice, one may use only the r > r c part of the axion profile obtained ignoring the core radius. Also, the typical momentum transfer in the monopole-nucleon scattering is much smaller than the mass of the vector boson, and therefore, the scattering cross section is sensitive to the axion profile far from the monopole.
What we want to know is the field configuration that minimizes the Hamiltonian. To this end, it is useful to change the radial coordinate r to a dimensionless variable z = r 0 /r. 5 Note that the theta term θF F in eq. (2.3) does not contribute the Hamiltonian density even if θ is spacetime dependent. The easiest way to prove this is to consider the following action where F µν = 1 2 µναβ F αβ / √ −g. Therefore, T µν θ = δS θ /δgµν = 0.

JHEP01(2020)185
The equation to solve is In order for the energy density of the axion-hidden monopole system to be finite, we also impose θ(z → ∞) = 0 and θ(z → 0) = θ 0 as the boundary conditions of eq. (3.4). Unfortunately, one cannot solve this differential equation analytically. Instead, let us solve this equation with its asymptotic forms as This has the analytical solutions given by Notice that this solution is smooth at z = √ m a r 0 , and we have checked that the asymptotic solutions given by (3.6) agree with the numerical solution of (3.4) to very high precision. See figure 1, where we make a comparison between the analytic and numerical solutions for θ 0 = 1 and √ m a r 0 = 10 −2 . One can see that the two lines agree very well. The monopole is located at z → ∞ (i.e. r = 0), and the axion field takes the vacuum value θ 0 at z → 0 (i.e. r → ∞). The axion field value starts to decrease around z 1 (i.e. r r 0 ), and this is due to the fact that a nonzero axion field costs large energy stored in the hidden electric field near the origin because of E H (r) 2 ∝ θ(r) 2 /r 4 . So far we have not specified the relationship between m a and 1/r 0 . For our purpose, we are interested in relatively light axion masses between 10 MeV and 10 GeV (the mass region of scalar particles searched by the beam-dump experiments), while the decay constant f H is constrained to be much larger than m a for this mass range. In the following calculation, we, therefore, assume m a r 0 1, for which the axion mass is relevant only in the tail of the axion field configuration around the monopole. As we shall see, those outer parts contribute most to the scattering cross section, because the typical momentum transfer between the monopole and the nucleon is smaller than m a . It is also clear that the precise size of the monopole core radius does not affect the calculation of the cross section.

Spin-dependent monopole-nucleon scattering cross section
With the axion field profile (3.6), we can now calculate the amplitude of the monopolenucleon scattering by introducing the following axion-nucleon interaction where N presents the nucleon (proton p or neutron n), C N is a constant of order unity, m N 1 GeV is the nucleon mass, and f a is a mass scale characterizing the strength of the axion-nucleon coupling. In the following calculation, we will assume f a = f H for simplicity, but f a can be much larger or smaller than f H in a slightly contrived setup. The values of C p and C n are actually model-dependent. In our numerical study, we choose |C p | 0.4 and |C n | 0.05, which are typical values in the DFSZ axion model [22].
In our set-up, the monopole is much heavier than nucleons, and so, we can treat a(x) as a classical scalar field around the monopole while we treat ψ N (x) as the usual quantized Dirac field in the same way as the Rutherford scattering. Then, the amplitude of the monopole-nucleon scattering mediated by the axion at the leading order is written as [23] iM M+N→M+N = C N m N θ(q) u N (p )γ 5 u N (p) , (3.8) where u N is the Dirac spinor, q = p − p is the transferred 3-momentum, and θ(q) is the 3-dimensional Fourier transform of θ(x) given by Here we have introduced the cut off in the integration corresponding to the core radius because the axion configuration derived before is valid only for r r c . Plugging eq. (3.6) into eq. (3.9), we have for √ m a r 0 > r 0 /r c . The integration of (3.10) and (3.11) involves an IR divergence at z → 0, which can be removed by inserting the regulator, lim δ→0 e −δ/z . Then, both integrations give approximately the same result, for √ m a r 0 1. 6 Here we have defined q = |q|/m a m N υ DM /m a 1, where υ DM 10 −3 is the typical velocity of DM with respect to the reference frame of the nucleon at the solar radius in our Galaxy. Also, we have numerically checked that eq. (3.12) gives a good approximation to eqs. (3.10) and (3.11) as long as √ m a r 0 1 no matter the size of r 0 /r c . This is because the DM-nucleon scattering is a low-energy process and therefore the integral is dominated by small z (i.e. large r).
For our benchmark point, the second term in the bracket of (3.12) is negligible since we consider very heavy monopole mass, namely the effect of the hidden monopole core can be neglected. Then, with eqs. (3.8) and (3.12), the resulting differential spin-dependent cross-section of the hidden monopole DM elastic scattering off a nucleon is evaluated as [23] dσ M+N→M+N dΩ We will use this result to constraint the axion decay constant by direct DM detection experiments in the next subsection. Before ending this subsection, let us briefly discuss the hidden vector boson DM in our model setup since it can also interact with nucleus through the axion exchange. The elastic scattering of the hidden vector boson DM and nucleon can be described by the following dimensional-7 effective operator as 14) and the resulting differential cross-section of the hidden gauge boson scattering off a nucleon is computed as Comparing eq. (3.15) with eq. (3.13), we find

16)
6 If ma 10 −3 GeV, then one has to do the integrations numerically, because q is no longer smaller than unity. For the interesting range of the decay constant, such a light axion is tightly constrained by astrophysics [14], and we do not consider it in this paper.

JHEP01(2020)185
and therefore, the hidden gauge boson-nucleon scattering is negligible. This motivates us to focus on the hidden monopole-nucleon scattering via the axion portal.
Lastly, with the results of (3.13) and (3.15) one may wonder why the hidden monopolenucleon scattering cross-section is suppressed only by the decay constant squared rather than fourth power of it. If one goes back to the calculation, one can see that the peculiar dependence on f a comes from the axion configuration around the monopole, and that the dominant contribution to the scattering cross section comes from r ∼ r 0 /m a , where the axion mass becomes relevant.

Limits and forecasts of direct DM search experiments
Now let us estimate the present limits on the predicted scattering cross section (3.13) from the direct DM search experiments. To this end, one can consider the following dimensional-6 pseudo-scalar operator given in ref. [24] since it gives the same dependence of the scattering cross section on the transferred momentum q, where χ is a Dirac fermionic DM with g χN being the coupling strength to the nucleon, and Λ is the cut-off scale of the theory. 7 The tree-level matrix element of the χ-N scattering in the non-relativistic limit is given by where m DM is the mass of χ, ξ(s) is the two-component spinor with s being the state of spin, and σ = (σ 1 , σ 2 , σ 3 ) is the Pauli spin matrix vector. The spin-averaged differential cross-section is then estimated as where we have assumed that m DM m N in the last equality. Thus, one can see that this has the same momentum dependence as (3.13). Note also that the operator (3.17) violates CP corresponding to the fact that the axion configuration around the monopole breaks CP.
In our study, we will focus on the case of proton since the COUPP and PICO experiments are sensitive to the DM scattering cross section with a proton. Making a comparison of (3.13) with (3.19) and referring the results analyzed in ref. [24], we present the constraints on the axion decay constant from the direct DM search experiments as a function of the axion mass with different choices of α H and m M in figure 4, where the amount of DM interacting with the nucleon is taken into account as explained in the previous subsection. In these plots, the red and green shaded regions are excluded by the current experiments COUPP-4 (4.0 kg CF 3 I) [25] and PICO-60 (52 kg C 3 F 8 ) [26], respectively, and the blue dashed line is the future prospect by PICO-500 (C 3 F 8 ) [27].

Implications for axion search experiments
In this section, we study the current limits and future prospects for the search of axions. In our setup (3.1), the axion is coupled not only to the SM particles but also to the hidden photons. The axion can decay into a pair of the hidden photons with the rate given by This additional decay mode changes various bounds on the axion decay constant from the beam-dump experiments (CHARM, SHiP, etc.) and some B-meson rare decays. We briefly JHEP01(2020)185 write down the relevant formulas and constraints in the following. For simplicity, we will set f H = f a . The beam-dump of CHARM (SHiP) experiment imposes that the number of events N det in the detector region is [32][33][34] where N a is the number of the axion produced in the solid angle of the detector with length ∆ = 35 (55) m, γ a = (1 − β 2 a ) −1/2 ≈ 10 (25) GeV/m a [34,35], c is the speed of light, and = 480 (70) m is the distance between the detector and the beam-dump. The lifetime of the axion τ a reads where Γ a is the total width of the axion, and Γ(a → γ H γ H ) and Γ(a → vis) are the rates of the axion decaying into the hidden photons and visible particles (e ± , µ ± , γ, · · · ), respectively.
On the other hand, the upper bound on the branching ratio of B 0 decaying into K 0 S plus missing energy is given by [32,36] where P esc (p a ) is the probability of the axion with momentum p a to escape from the detector with B(a → γ H γ H ) = Γ(a → γ H γ H )/Γ a and B(a → vis) = Γ(a → vis)/Γ a being the branching ratios of the axion decaying into the hidden photons and visible particles, respectively. In our numerical study, we take max = 4 m as required in ref. [32]. Applying eq. (4.1) to (4.5), we present in figure 3, how the decay channel a → γ H γ H affects the limits from the CHARM experiment, the SHiP experiment, and B 0 decaying into K 0 S with missing energy. The constraints from B → Kµ + µ − and B → µ + µ − are not shown here since we have checked that the effects of them are weaker than the others. Notice that we assume the Yukawa-like coupling for the axion in our analysis, where the couplings between the axion and the SM fermions are proportional to the SM Yukawa couplings [32]. This interaction may be induced in a set-up like the DFSZ axion model [37,38].
In the upper plots of figure 3, the parameter region of f −1 a above the color lines corresponds to the case where the axion decays before arriving the detector. Due to the JHEP01(2020)185 hidden photon decay mode of the axion, the upper parts of these color lines shrink as one has to reduce the coupling strength (∝ f −1 a ) between the axion and the SM particles in order to maintain the same predicted event numbers. On the other hand, the parameter region of f −1 a below the color lines corresponds to the case where the axion decays after leaving the detector. In this region, the lifetime of the axion is very long such that the predicted event numbers is N a ∆ X=e,µ,γ Γ(a → XX)/(γ a β a c), which is independent of α H . This explains why the lower parts of the color lines are insensitive to α H . In the bottom plot of figure 3, the excluding region of f −1 a is enlarged since we identify with the axion decays into the hidden photons within the detector as an escaped event, then B(a → γ H γ H ) gives a non-negligible contribution in eq. (4.5) even it is a sub-dominant decay mode at m a > 2m µ .
Combining The latter is opened compared to the previous analysis [32], and this is due to the hidden photon decay channel of the axion. Similar results would be obtained with other benchmark points.
Lastly, let us comment on the SN1987A supernova bound. The QCD axions (and other light particles with feeble interactions) are considered to be produced in the proto-neutron star (see refs. [39,40] and references therein). If the QCD axion carries away too much energy from it, the neutrino burst duration would be shorter than ∼ 10 sec. in contradiction with the observation. 8 Compared to the QCD axion, the axion in our scenario has a much heavier mass, and it decays into hidden photons with a sizable branching fraction. The temperature of the proto-neutron star can be as high as a few tens MeV, and so, the parameter region (i) may be subject to the SN1987A bound. 9 A more detailed analysis is necessary to estimate the bound. In any case, the SN1987A supernova bound is not included in our analysis partly due to the incomplete understanding of the supernova explosion mechanism.

Discussion and conclusions
First, let us estimate the hidden photon contribution to dark radiation. In the minimal case without the axion portal coupling, the number of the extra neutrino species, ∆N eff , was estimated in ref. [8], where the entropy stored in the W ± H is transferred to the hidden photons after the temperature falls below the mass of W ± H . Let us see how this is modified in our model setup. To be concrete, we fix the axion mass and decay constant to be the region of   keeps the hidden photons in equilibrium with the SM sector at temperatures well below the mass of W ± H or the PQ fermions, and the hidden photon is considered to decouple from the SM plasma at a temperature around the electroweak scale. Therefore, we do not need to take account of the entropy stored in the heavy hidden particles. For the reference parameter, the axion remains in equilibrium and abundant until the temperature becomes comparable to the axion mass. In particular, the axion mainly decays into the hidden photons when the axion becomes non-relativistic and its abundance gets Boltzmann suppressed. This production of the hidden photons dominates over the contribution at the decoupling of the hidden photons. Hence, we have ∆N eff 4/7 + 0.05 ∼ 0.6, which is consistent with the estimation given in ref. [42] assuming the presence of the axion coupling JHEP01(2020)185  to leptons. In the other viable parameter region, m a = O(100) MeV, f a = O(10 4 ) GeV, the contribution to dark radiation will be smaller because the axion disappears from the plasma at higher temperatures. The dark radiation with ∆N eff 0.5 can relax the H 0 tension significantly [15]. 10 Next let us comment on a kinetic mixing between U(1) H and the hypercharge U(1) Y . As discussed in refs. [43,44], we can write the following dimensional-5 operator where B µν is the field strength tensor of U(1) Y , and M is a cut off scale. This effective operator gives rise to the kinetic mixing of υ H /M , which is of order 10 −13 for υ H = 10 5 GeV and M = 10 18 GeV. Then, both monopoles and W ± H acquire a fractional electric charge of order υ H /M through the kinetic mixing. The current bound on such mini-charged DM is satisfied for the mass of O(100) TeV and the kinetic mixing of O(10 −13 ) [46,47].
In this paper, we have studied the possibility that the hidden monopoles which account for a sizable fraction of the observed DM density are coupled to the SM sector via the axion portal coupling. We have determined the axion configuration around the monopole so that JHEP01(2020)185 it minimizes the energy of the system, and then estimated the scattering cross section of the monopole with a nucleon via the axion portal. Using the fact that the dependence of the cross section on the transferred momentum is the same as that for a CP-violating pseudo scalar coupling between quarks and a Dirac DM particle, we have translated the experimental bounds on the spin-dependent scattering cross section with a nucleon to limit the model parameters. The axion itself can also be searched for at the beam-dump experiment. Combining the current limits and future prospects for the search of DM and axions, we have found two viable regions around