Probing pseudo Nambu-Goldstone boson dark matter at loop level

In the standard model extended by a single complex scalar boson with a softly broken global U(1) symmetry, a pseudo Nambu-Goldstone boson becomes a candidate for dark matter. In this paper, we discuss the direct detection of the pseudo Nambu-Goldstone boson dark matter. Since the tree-level amplitude for dark matter-nucleon scattering vanishes, higher order quantum corrections for the amplitude should be taken into account. We perform the calculation at the next-to-leading order in QCD in a systematic manner.


Introduction
Although it is well-known that non-baryonic dark matter (DM) exists in the universe, nature of dark matter is still a mystery of the universe except the fact that its relic abundance occupies about 26% of the total energy density of the universe [1,2]. A prominent candidate for DM is a stable and non-relativistic particle that has the weak scale interactions. In this case, the annihilation rate determining its relic abundance is closely correlated with scattering rates with the Standard Model (SM) particles and production rates at collider experiments through the crossing symmetry.
Direct detection experiments of cold dark matter explores scattering events between DM and nuclei. So far, any viable signal of DM has not been found even in recent ton-scale detector experiments, which leads to the bounds on the interaction between DM and nucleon. The most stringent upper bound on the spin independent (SI) elastic scattering cross section between DM and nucleon is given by the XENON1T Collaboration [3], which is, for example, 4.1 × 10 −47 cm 2 at the DM mass of 30 GeV at 90% confidence level.
The current null results of direct signal of DM motivate us to consider a framework where the interactions between DM and nucleon are suppressed in non-relativistic limit. One of the ideas is to consider a pseudo Nambu-Goldstone DM [4] or a pseudo scalar portal fermionic DM [5,6,7]. In the former case all the interactions between DM and SM particles are described by derivative couplings at the tree level, while in the latter case the DM spinor product for elastic scattering is proportional to the DM velocity in non-relativistic limit. As a result, the elastic scattering amplitude is necessarily suppressed by the small DM velocity at tree level, thus these DM candidates can naturally be consistent with the strong constraints from the current direct detection experiments.
In this paper we consider a pseudo Nambu-Goldstone DM proposed in ref. [4], and study the possibility to detect it directly. As stated above, the elastic scattering amplitude for direct detection is suppressed at the tree level in this model, and it can be exactly zero in non-relativistic limit. However, the scattering amplitude is expected to be finite at the loop level. We will perform the calculation at one-loop level for non-QCD part and at the next-to-leading order level in QCD based on the formalism given in ref. [8] where the scattering processes with gluon in nucleon, which are sometimes missed in the literature, are systematically taken into account.
Although the next-to-leading order calculation in QCD is not necessary for a rough evaluation, it gives O(10%) corrections in the amplitude and the theoretical uncertainty regarding perturbative QCD calculation is reduced significantly. We explore a parameter space consistent with the observed DM relic abundance, the SM Higgs boson decay, and the perturbative unitarity bound. We also compare the predicted elastic scattering cross section with the sensitivity of the future direct detection experiment DARWIN [9].

The Model
We consider the SM augmented by a complex scalar field S with a softly broken global U (1) symmetry. The model is invariant under the transformation S → e iα S with a real constant α except for the soft breaking term. The scalar potential of the model is given by where H is the SU (2) L Higgs doublet which couples to the SM particles. The last term corresponds to the soft breaking term of the global U (1) symmetry. Due to the tachyonic mass terms for the scalar fields, both H and S acquire vacuum expectation values (VEVs) in a wide range of the parameter space, which is the situation we are interested in. Then these fields are expanded around the vacuum as where G + and G 0 are the Nambu-Goldstone bosons associated with the electroweak symmetry breaking, v ( 246 GeV) and v s are the VEVs for H and S, respectively.
h and s are the CP even scalar fields while χ is the CP odd scalar field which is the would-be Nambu-Goldstone boson. Due to the soft breaking term of the global U (1) symmetry, non-zero mass for χ arises. Even after the symmetry breaking, a Z 2 symmetry remains, which stabilizes χ and makes it a candidate for DM.
Due to the symmetry breaking, the CP even states h and s mix via the Higgs portal coupling λ HS . We will derive the mass eigenstates h 1 and h 2 at one-loop level.
As mentioned in Introduction, χ-q scattering amplitude vanishes at the tree level.
Therefore, one-loop corrections are necessary for studying direct detection of χ DM.
To this end, we execute the calculation following ref. [10]. In the literature, the inverse propagators for the scalars are calculated diagramatically. Then the mass matrices for the CP odd and even sectors are defined by taking zero external momenta. That corresponds to the one obtained from the effective potential. Then the mass eigenstates are given by diagonalizing the mass matrices. In the following calculation we adopt Landau gauge and MS renormalization scheme as in the literature. All couplings, scalar fields, and VEVs are renormalized values. We will see that the renormalization scale dependence is cancelled in the amplitude as expected.
The mass matrices for the CP even and CP odd sectors are given by where T h and T s are renormalized tadpoles for h and s, which satisfy the stationary conditions: Then, the inverse propagators are given as where indices i, j represent h, s (or h 1 , h 2 ), G 0 , χ, and Π ij (p 2 ) correspond to the renormalized self-energies with the external lines i and j. The concrete expressions for Π ij (p 2 ) are collected in Appendix. Here we have introduced the quantities ∆Π ij (p 2 ) defined by 8) and the definition ofM 2 ij in Eq. (2.7) follows accordingly. Note thatM 2 ij correspond to the mass matrix derived from the effective potential, i.e., zero external momenta.
Since the mass matrix M 2 odd is diagonal, the physical (pole) masses for the CP odd fields are simply given by Γ G 0 G 0 (0) = 0 and Γ χχ (m 2 χ ) = 0 where m χ is the pole mass of χ, i.e., The CP even sector, on the other hand, needs to be diagonalized. Following ref. [10], we derive an eigenstate basis with one-loop correction by diagonalizingM 2 ij . Using the equations given above,M 2 ij in the CP even sector are rewritten as 14) Then the mixing angle for the diagonalization is obtained by The eigenstates h 1 h 2 are then given by (h 1 , h 2 ) T = O T (h, s) T . We define h 1 as the lighter field, which is identified as the observed Higgs boson. Their physical masses are then derived straightforwardly as Using the mass eigenstates, the scalar potential can be expanded around the VEVs. In the mass eigenstate basis, the coefficients for scalar cubic and quartic couplings, which are relevant for our discussion, are expressed as where i, j, k, m are 1 or 2. Additionally, the Yukawa couplings to quarks are given by where q are quarks.

Scattering Cross Section
In this section we compute the SI cross section of χ DM with nucleon. To avoid confusing readers, we clarify some terminology related to our calculation. We will perform the calculation literally at one-loop level for non-QCD related part. For QCD part, on the other hand, the amplitude is derived at the next-to-leading order (NLO) level. Throughout this paper we use the term NLO or leading order (LO) in terms of order of QCD strong coupling α s . For example, LO contains oneloop diagrams for χ-g scattering. However, the gluon contributions to the effective scalar coupling are O(α 0 s ) [8,11,12]. That is why we use the term LO for such χ-g processes, and similar discussion is applied for NLO. As we will see, the gluon contributions become important in some parameter space.

Formalism
We briefly summarize the formalism for the calculation of the SI scattering cross section of a real scalar DM with nucleon based on refs. [8,13]. Using the formalism in ref. [8], we calculate the scattering amplitude at the NLO in QCD.
The effective Lagrangian relevant for the scattering process is Here G a µν represents the field strength tensor of gluon field and the quark masses are denoted as m q . The operators O q µν and O G µν are the twist-2 operators of quarks and gluon, respectively, which are defined by with D µ the covariant derivative. Then the SI scattering cross section of χ with nucleon N is obtained as [8,13] where m N is the nucleon mass, and f N scalar and f N twist2 are given by Here f N T q , f N T g , q N (2; m Z ),q N (2; m Z ), and g N (2; m Z ) are the matrix elements of the effective operators in nucleon state. µ had is the hadronic scale (i.e., around 1 GeV), and m Z is the Z boson mass. The numerical values for these quantities are given in ref. [8] #1 based on the QCD lattice simulation [14,15] and CETEQ-Jefferson Lab collaboration [16]. As we will see, the contribution to the twist-2 type operators is negligibly small. Therefore, the SI cross section is determined by the scalar-type interactions.

Wilson coefficients
First of all, we derive the effective Lagrangian from full theory by matching at the weak scale denoted as µ W m Z . As described in Introduction, the tree-level amplitudes for χq → χq are cancelled in the non-relativistic limit. Therefore loop-level calculations are necessary to evaluate the scattering amplitude in the limit. There T q in ref. [8] and so on. We have additionally introduced f N T g defined as (−8/9)f N T g ≡ N | αs π G a µν G aµν |N /m N evaluated at three flavors, which leads to Figure 1: Feynman diagrams for χ-q and χ-g scattering processes.
are three types of diagrams shown in Fig. 1; (i) Self-energies, (ii) Vertex corrections, and (iii) Box and triangle diagrams. #2 The most part of the computation for the diagrams (i) has already been done in the previous section. Let us apply the results to the scattering process. The diagrams (i) give rise to scalar-type interactions. The matching at the weak scale for q = u, d, s, c, b where ∆ 12 is given in Eq. (2.15), and∆ 22 is given bỹ The computation of the diagrams (ii) is rather simple. These diagrams give vertex corrections to the χ-χ-h i couplings. Denoting them as ∆c χχh i (collected in Appendix), the Wilson coefficients are obtained as

11)
#2 The NLO diagrams for QCD part are not depicted for simplicity, but these are taken into account in the numerical study.
for q = u, d, s, c, b. It is noted that, although both C q S (µ W )| self and C q S (µ W )| vert have the renormalization scale dependence, the dependence is cancelled in C q S (µ W )| self + C q S (µ W )| vert as expected.
The χ-q scattering process drawn in the diagrams (iii) gives both scalar and twist-2 type contributions. However, the resultant Wilson coefficients are proportional to , d, s, c, b), thus they are negligibly small. For the χ-g scattering, on the other hand, the top loop diagram should be taken into account since there is no such suppression. They can be calculated easily in Fock-Schwinger gauge (see Appendix for details) and the resultant expressions are We refer to the terms proportional to J ij box and J ij tri as 'box' type and 'triangle' type, respectively. Note that there is no m χ dependence in J ij tri /m 2 χ , which is obvious from the corresponding diagrams. It is found numerically that the contribution to the amplitude from the box-type diagrams is much smaller than that from the triangletype diagrams in the parameter space we are interested in. Here we have ignored NLO contributions and we will treat it as a theoretical uncertainty as in ref. [8]. This is because it is expected to be suppressed compared to the other NLO contributions.
To summarize, the weak scale matching gives (3.14) The Wilson coefficients at the hadronic scale µ had are obtained by the renormalization group equations, along with the matching at bottom and charm mass scales, consistently at the NLO in QCD [8].

Numerical results
Now we are ready to show the numerical results. The SI cross section of χ DM with proton is plotted in Fig. 2 [18]. The grey region is excluded by the perturbative unitarity bound λ S ≤ 8π/3 [19]. The red band represents the parameter space which can reproduce the observed DM relic abundance within 3σ range of the PLANCK Collaboration data [2] in the thermal freeze-out scenario. It has been found that since the SI cross section is so suppressed in a wide range of parameter space that there is no substantial constraint on the model from the current experimental limit of the XENON1T experiment [3]. To be concrete, the region excluded by XENON1T is always in the unitarity bound. In the plots, the future reach of the DARWIN experiment [9] is also shown in the orange dot-dashed line, assuming that all the DM abundance is composed by χ. It indicates that a part of the parameter space (70 GeV m χ 100 GeV) can be probed by the DARWIN experiment, where the thermal relic DM scenario can reproduce all the DM abundance. We have seen a similar indication for larger sin θ values. Therefore the future direct detection experiments will be able to probe a part of the thermal freeze-out scenario with m χ ∼ 100 GeV.

Conclusion
We have studied the detectability of DM in a model where a complex singlet scalar is added to the SM. In the model, a softly broken global U (1) symmetry has been assumed, and a would-be Nambu-Goldstone boson χ becomes a candidate for DM due to a remnant Z 2 symmetry. Since χ interacts with the SM particles via socalled Higgs portal coupling, it would be possible to directly detect χ DM via χnucleon scattering. It is known, however, that the tree-level scattering amplitude vanishes in non-relativistic limit. Thus we have taken into account the scattering process at one-loop level for non-QCD part. For QCD effect, on the other hand, the scattering amplitude has been calculated at the next-to-leading order in QCD strong coupling systematically to reduce the theoretical uncertainties. It has been found that the predicted SI cross section is small in a wide range of the parameter space of the model, and there is no substantial bound from the current direct detection experiments. However, a part of the parameter space, which includes canonical thermal relic scenario accounting for the present DM abundance, will be probed in the future direct detection experiment DARWIN. Note added: While completing our paper, we found that ref. [21] studied the direct detection of dark matter in the same model. We agree qualitatively with their results in large dark matter mass region, as well as in the other dark matter mass region if the diagrams (iii) are ignored. On the other hand, we have seen a different behavior in low dark matter mass region. This is due to the diagrams (iii), which are discarded in their study. Although these diagrams are irrelevant for the estimation of the current bound from the XENON1T experiment, it will be important for the future study of this model in the direct detection experiments.

A Loop functions
The loop functions are basically expressed by so-called A 0 function and B 0 function (and their derivatives), .
where P 2 = m 2 χ in the last two functions. Here the divergent pieces are subtracted in the MS renormalization scheme implicitly. We use the analytic expressions for the loop functions, which are numerically checked by using LoopTools [20].

A.1 Tadpoles
The tadpoles T h i are given by A 0 function as, Then T h and T s are obtained by rotating with the orthogonal matrix O given in

A.3 Vertex corrections
The corrections to the cubic couplings c χχh i are given by As described in the main text of the paper, the renormalization scale dependence is cancelled in the Wilson coefficients, C q S (µ W )| self + C q S (µ W )| vert . It would be helpful to see how it behaves in small and large m χ limit. For m χ → 0 limit, we have found that it is proportional to m 2 χ . For example, #4 In large m χ limit, on the other hand, its absolute value increases logarithmically, e.g.,

A.4 Box and triangle diagrams
To compute the Wilson coefficients induced by χ-g scattering, it is legitimate to derive Higgs correlation functions shown in Fig. 4. We denote them asΠ ij (q 2 ), which are obtained straightforwardly by using the formula given in refs. [11,12] as, where c G (q 2 ; µ W ) = − c qqh i c qqh j 8m 2 χ I(−q 2 /m 2 χ , x t ) . (A.33) Here I(t, x t ) is a dimensionless function, with β(t, x t ) = 1 + 4x t /t, x i = m 2 h i /m 2 χ , and x t = m 2 t /m 2 χ . Then, the Wilson coefficients coming from the box and triangle diagrams are given by the following integrals, with κ ij = 2 for i = j otherwise 1. As mentioned in the main text of the paper, J ij tri /m 2 χ is constant with respect to m χ . On the other hand, J ij box ∝ m 4 χ (m 3 χ ) for small (large) m χ region. Therefore the contribution from the box diagrams is suppressed as 1/m χ in large DM mass region while it becomes constant in small DM mass region.