# Theory of record thermopower near a finite temperature magnetic phase
transition: IrMn

Peter Wölfle  Institute for Theory of Condensed Matter, Karlsruhe Institute of
Technology, 76128 Karlsruhe, Germany  Institute for QuantumMaterials and
Technologies, Karlsruhe Institute of Technology, 76021 Karlsruhe, Germany
Timothy Ziman  Institute Laue-Langevin, 71 Avenue des Martyrs, 38042 Grenoble,
France

###### Abstract

The effect of scattering of conduction electrons by dynamical spin
fluctuations on the thermopower in metals near a thermal phase transition into
an antiferromagnetic phase is considered. We are interested in a transition at
room temperature, as has been studied in a heterostructure involving layers of
IrMn. We show that the electrical resistivity exhibits a narrow but low peak
at the transition, which may be ddifficult to detect on top of the main
contributions induced by phonons and impurities. By contrast, the thermopower
is found to exhibit a prominent peak both as a function of temperature $T$ for
fixed layer thickness $t_{AFM}$ and as a function of $t_{AFM}$ for fixed $T.$
We conjecture that the transition temperature $T_{c}$ is a function of both
$t_{AFM}$ and the Fermi energy $\epsilon_{F}$. Both dependencies give rise to
a sharp peak of the thermopower as a function of $T$ or $t_{AFM}$ near the
transition. The estimated magnitude of the peak for the case of three-
dimensional longitudinal spin fluctuations is in good agreement with
experiment.

##  I Introduction

The recently observed ftacular temperature dependence of the thermopower in a
magnetic heterostructure involving the antiferromagnetic metal IrMn Ziman20
suggests that the spin fluctuations near the antiferromagnetic transition may
be responsible for the observed peak at the ordering temperature. This finding
is all the more interesting in that the magnetic transition temperature may be
tuned by the thickness of the IrMn layers to be at room temperature, which
makes the effect highly promising for applications. The effect of elastic
scattering of conduction electrons by the local spins of magnetic metals near
the phase transition into a magnetically ordered phase has been studied first
by DeGennes and Friedel DeGennes58 , who assumed the spin configuration to be
temperature dependent and static as given by an Ornstein-Zernike form.
Somewhat later, Fisher and Langer Fisher68  revisited the problem in the light
of the theory of classical critical phenomena. These authors pointed out that
the equal time spin correlation function entering the scattering cross section
also appears in the internal magnetic energy and is therefore related to the
specific heat. They also mentioned that inelastic scattering may be important.
The effect of inelastic scattering on the thermopower was considered by Entin-
Wohlman, Deutscher, and Orbach Entin76  in a model calculation leaving the
dynamics of the local spin system as an input quantity to be determined from
case to case. A detailed model calculation of the electrical resistivity of
antiferromagnetic metals, in the framework of the Self-Consistent
Renormalization Theory of spin fluctuations in itinerant magnets Hasegawa74 ;
Moriya85  has been worked out by Ueda Ueda77 . In all these previous studies
the anomalies near the transition found for the transport properties appeared
to be relatively weak and cannot account for the prominent peak found in the
thermopower as a function of temperature or as a function of layer thickness
of an IrMn heterostructure Ziman20 .

In this paper we estimate the contribution of scattering of the charge
carriers by dynamical spin fluctuations by assuming a phenomenological form of
the spin excitation spectrum dictated by symmetry and conservation laws. We
observe that the presence of gapless fermionic excitations in a metal changes
the spin excitation spectrum of the local spins in a decisive way, leading to
a strongly temperature dependent, at $T_{c}$ seemingly divergent contribution
to the resistivity. The divergent behavior arises for not too high transition
temperature, $T_{c}\ll\epsilon_{F}$, where $\epsilon_{F}$ is the Fermi energy
(here and in the following we use energy units of Kelvin). The divergence is
cut off close to $T_{c}$ at a temperature $T_{x}$, marking the transition into
a quantum critical regime. The effect may be characterized as a continuation
of quantum critical scattering in a narrow region of the phase diagram along
the phase boundary. There are two effects introduced by the itinerant
electrons into the spin dynamics of antiferromagnetic metals as contrasted to
insulators. First, the transition temperature is shifted by an amount
proportional to the static wave vector dependent conduction electron spin
susceptibility at the ordering wave vector $\mathbf{Q}$,
$\chi_{s}(\mathbf{Q,}0)$. The latter depends on the Fermi energy and thus
gives rise to a dominant contribution to the thermopower. Secondly, the
dynamics of the spin fluctuations at low energy is dominated by the Landau
damping mechanism leading to inelastic scattering processes strongly enhanced
near the transition.

##  II Model and Method

###  II.1 Hamiltonian

We assume a system of interacting localized spins, coupled to conduction
electrons as expressed by the Hamiltonian

$H=H_{c}+H_{S}+H_{ex}$ (1)

where $H_{c}$ represents a single conduction band

$H_{c}=\sum_{\mathbf{k},\alpha}\epsilon_{\mathbf{k}}c_{\mathbf{k}\alpha}^{{\dagger}}c_{\mathbf{k}\alpha},$
(2)

We assume the conduction electron system to be three-dimensional.

The dynamics of the localized spin system is defined by

$H_{S}=\sum_{\mathbf{q}}\sum_{\alpha=x,y,z}I_{\alpha}S_{\mathbf{q}}^{\alpha}S_{-\mathbf{q}}^{\alpha},$
(3)

where we allow for anisotropic interaction $I_{z}\gg I_{x,y}$ caused by the
strong spin-orbit interaction at the Ir ions.

The coupling of conduction electrons to localized spins is described by

$H_{c-S}=J\sum_{\mathbf{q}}\mathbf{s}_{\mathbf{q}}\mathbf{S}_{-\mathbf{q}}.$
(4)

Here $S_{\mathbf{q}}^{\alpha}$ are the Fourier components of the localized
spin operators, and the conduction electron spin operator is defined by
$\mathbf{s}_{\mathbf{q}}=\sum_{\mathbf{k},\alpha,\beta}\mathbf{\tau}_{\alpha\beta}c_{\mathbf{k+q}\alpha}^{{\dagger}}c_{\mathbf{k}\beta}$
, with $\mathbf{\tau}_{\alpha\beta}$ the vector of Pauli matrices (the
coupling constants $I_{\alpha},J$ are given in units of [energy/density]).

###  II.2 Spin fluctuations in the paramagnetic phase

The transport properties of the conduction electron system of MnIr in the
temperature range around room temperature are governed by electron-phonon
interaction and the exchange interaction mediated by $H_{c-S}$ as we argue
below. We assume the system to be anisotropic in spin space, with preferred
direction along the $z-$axis. We will consider both three-dimensional and two-
dimensional spin fluctuations, where the $3d$ model appears to describe the
experiment Ziman20  better, as we shall see.

In the absence of coupling of local spins and conduction electron spins, for
$J=0$, the longitudinal susceptibility of localized spins $\chi_{loc}^{zz}$ is
assumed to be well approximated by the static Ornstein-Zernike form

$\chi_{loc}^{zz(0)}(\mathbf{q,}\omega)\approx\frac{n_{S}}{I_{z}}\frac{1}{(T-T_{c,I})/T_{c,I}+(\mathbf{q-Q})^{2}\xi_{0}^{2}}$
(5)

where $\mathbf{Q}$ is the ordering wave vector and $\xi_{0}$ is a microscopic
spin interaction length of the order of a lattice spacing and $n_{S}$ is the
density of localized spins. The transition to the antiferromagnetic phase is
signaled by the divergence of $\chi_{loc}^{zz(0)}$ at wave vector
$\mathbf{q=Q}$ and at temperature $T_{c,I}=O(I_{z})$. In the quasi two-
dimensional slab geometry of IrMn within the nanostructure studied in Ziman20
the transition temperature is found to depend on the layer thickness
$t_{AFM}$. One source of such a dependence is the quantization of the momentum
component $(\mathbf{q-Q})_{z}$ normal to the layer surface,
$(\mathbf{q-Q})_{z}=\frac{\pi}{t_{AFM}}(2n+1)$, $n=0,\pm 1,...$. The minimal
value of $(\mathbf{q-Q})_{z}$ at $n=0$ gives rise to a suppression of the
transition temperature

$T_{c,I}^{\ast}=T_{c,I}[1-(\frac{\pi\xi_{0}}{t_{AFM}})^{2}]$ (6)

Figure 1: Transtion temperature into the antiferromagnetic phase versus layer
thickness $t_{AFM}$. Theory according to Eq. (6) (solid line). Data points
from Ziman20 : Seebeck measurements (red dots); susceptibility data (blue
dots); XMLD data (green dots).

In Fig. 1 we show that the $T_{c}$ data of Fig. 4 of Ziman20  may be fitted
reasonably well by Eq. (6) using a spin interaction length $\xi_{0}=0.17$nm
and $T_{c,I}^{\ast}=308$K. Here the reference transition temperature
$T_{c,I}^{\ast}$ is not the bulk value, which is much higher, but a reduced
value appropriate for the composite layer structure. For example, the
effective spin exchange energy (the quantity $I_{z}$ in our notation) has been
found to be strongly varying with the Fermi energy Umetsu04 , possibly
accounting for a reduction of the bulk $T_{c}$ value by a factor of $\sim 3$.
The above estimate differs somewhat from the theoretical results for $T_{c}$
versus $t_{AFM}$ given in Frangou16  and the model calculation for the
suppression of $T_{c}$ based on the reduction of the interaction volume of a
spin near the surface Zhang01 .

The conduction electron spin susceptibility in the model of noninteracting
fermions is given at low energy by

$\chi_{s}(\mathbf{q,}\omega)\approx\chi_{s}(\mathbf{q,}0)[1+i\pi\frac{\omega}{v_{F}q}]$
(7)

where $v_{F}$ is the Fermi velocity (we employ an isotropic band structure for
simplicity) and we assume isotropy in spin space. In the presence of coupling
the susceptibility of localized spins is renormalized approximately as

$\chi_{loc}^{zz}(\mathbf{q,}\omega)\approx\frac{1}{[\chi_{loc}^{zz(0)}(\mathbf{q,}\omega)]^{-1}+J^{2}\chi_{s}(\mathbf{q,}\omega)}$
(8)

In the vicinity of the transition and for small $\omega$ the susceptibility
$\chi_{loc}^{zz}$ is sharply peaked at $\mathbf{q\approx Q}$ so that
$\chi_{s}(\mathbf{q,}\omega)$ may be replaced by
$\chi_{s}(\mathbf{Q,}\omega)$, and $\chi_{loc}^{zz}$ takes the form

$\chi_{loc}^{zz}(\mathbf{q,}\omega)\approx\frac{n_{S}}{I_{z}}\frac{1}{\tau+(\mathbf{q-Q})^{2}\xi_{0}^{2}+\frac{i\omega}{\gamma}}$
(9)

Here we defined the reduced temperature $\tau=(T/T_{c}-1)$ and a renormalized
transition temperature $T_{c}=T_{c,I}^{\ast}-\Delta T_{c}$ where $\Delta
T_{c}=T_{c,I}^{\ast}n_{S}J^{2}\chi_{s}(\mathbf{Q,}0)/I_{z}\approx
T_{c,I}^{\ast}(n_{S}J^{2}/I_{z})N_{F}$, with $N_{F}=mk_{F}/\pi^{2}$ the
conduction electron density of states at the Fermi level ($m$ is the mass and
$k_{F}$ is the Fermi wave number of the conduction electrons). The transition
temperature $T_{c}$ is seen to depend on the Fermi energy
$\epsilon_{F}=k_{F}^{2}/2m$ , a fact of considerable importance for the
thermopower as we show later. The reference energy $\gamma$ entering the
Landau damping term in the denominator of $\chi_{loc}^{zz}(\mathbf{q,}\omega)$
is defined as
$\gamma=[(n_{S}J^{2}/I_{z})\chi_{s}(\mathbf{Q,}0)]^{-1}v_{F}Q/\pi\approx(I_{z}/\pi
n_{S}J^{2}N_{F})v_{F}Q\approx T_{c}(N_{F}J)^{-2}$, assuming $Q\approx k_{F}$
and $n_{S}\approx n\approx N_{F}\epsilon_{F}.$

###  II.3 Spin fluctuations in the antiferromagnetic phase

In the ordered phase the translation invariance is broken as indicated by the
ordering vector $\mathbf{Q}$. The rotation invariance in spin space is also
broken by the appearance of the ordered moment vector $\mathbf{M(R}_{j})$ .
Spin fluctuations along $\mathbf{M}$ (longitudinal) and orthogonal to
$\mathbf{M}$ (transverse) behave in a fundamentally different way. The
magnetic order also affects the quasiparticle excitations in major ways. The
Brillouin zone gets shrunk into the magnetic Brillouin zone, leading to
multiple bands induced by back-folding. To keep things simple we will assume
two-sublattice antiferromagnetic order ($\mathbf{Q=(}\pi,\pi,\pi)$) for which
case one finds two bands, dubbed ”valence” and ”conduction” band, separated by
a gap $\Delta$. The structure of the quasiparticle and collective excitations
has been calculated for the Hubbard model within the Random Phase
Approximation Moriya12 ; Schrieffer89 ; Chubukov92 .

####  II.3.1 Longitudinal spin fluctuations

The longitudinal spin susceptibility is finite in the ordered phase but
diverges upon approach to the thermal transition. This is expressed by the
form

$\chi_{loc}^{zz}(\mathbf{q},\omega)\approx\frac{n_{S}}{I_{z}}\frac{1}{r_{l}+(\mathbf{q-Q})^{2}\xi_{0}^{2}+i(\omega/\gamma)},$
(10)

where $r_{l}=c_{l}M_{s}^{2}=c_{l}M_{0}^{2}|\tau|$, assuming mean field
behavior of the ordered moments $M_{s}\propto M_{0}|\tau|^{1/2}$, where
$M_{0}$ is the saturated moment. Here $c_{l}$ is a constant of $O(1)$.

####  II.3.2 Transverse spin fluctuations

If spin rotation invariance in the plane orthogonal to the ordered moments
$\mathbf{M(R}_{j})$ is still preserved, which we suspect not to be the case
for MnIr, gapless spin excitations (spin waves) exist. In this case the
transverse spin susceptibility involving spin excitations orthogonal to
$\mathbf{M(R}_{j})$ is divergent in the whole ordered phase, not just at the
critical point. More generally, the transverse spin susceptibility may be
expressed as

$\chi^{-+}(\mathbf{q},\omega)\approx\frac{n_{S}}{I_{z}}\frac{1}{\Delta_{s}/\gamma+(\mathbf{q-Q})^{2}+i(\omega/\gamma)}$
(11)

where $\Delta_{s}$ is the spin wave gap. The coupling of the transverse spin
fluctuations to the quasiparticles depends on the respective quasiparticle
bands. For definiteness we assume that the Fermi energy lies in the lower
(valence) band. Then the scattering of quasiparticles from the valence band
into the valence band is suppressed, the corresponding vertex
$\Lambda_{vv}^{tr}(\mathbf{k,k+q+Q)}\propto\mathbf{q\cdot\nabla}\epsilon_{\mathbf{k}}$
such that the product $(\Lambda_{vv}^{tr})^{2}\chi^{-+}$ appearing in the
expression for the self energy is no longer singular and hence does not give
rise to critical behavior Schrieffer89 ; Chubukov92 ; Adler65 . By contrast,
the vertex function for scattering from the valence band into the conduction
band is of $O(1)$. However, the latter excitation requires a minimum energy
$\Delta$ and is thus possible only in the vicinity of the transition when
$T\gtrsim\Delta$. In the following we will assume that the spin wave gap
$\Delta_{s}>T_{c}$ that spin waves may not be thermally excited.

##  III Electrical resistivity

We now present a model calculation of the contribution of scattering by
antiferromagnetic fluctuations to the electrical resistivity $\varrho$. Since
the typical momentum transfer in such scattering processes is large, of order
of the ordering wave vector $\mathbf{Q}$, and therefore, assuming a half-
filled conduction band, of order $k_{F}$, the momentum relaxation rate
$1/\tau_{tr}$ is approximately related to the imaginary part of the electron
self energy $\Sigma(\mathbf{k,}\omega)$ by $1/\tau_{tr}\approx
2\mathit{Im}\Sigma$. The resistivity is then given in terms of the self energy
as

$\rho\approx\frac{m}{e^{2}n}2\mathit{Im}\Sigma(k_{F}\mathbf{,}\omega\approx
T)$ (12)

###  III.1 Three-dimensional spin fluctuations

In a clean metal at low temperature scattering by antiferromagnetic
fluctuations affects only a small part of the Fermi surface: the ”hot spots”
connected by the ordering wave vector $\pm\mathbf{Q}$. At higher temperatures
phonon scattering, or at all temperatures impurity scattering, helps to remove
the constraints imposed by momentum conservation so that the critical behavior
at the hot spots is distributed all over the Fermi surface. The self-energy
may be calculated to one-loop order, taking into account phonon and impurity
scattering, which leads to a prefactor $A=O(1)$.

In the case of electrons and spin fluctuations that are both three-
dimensional, we get

$\displaystyle\mathit{Im}\Sigma(\mathbf{k,}\omega)$ $\displaystyle\sim$
$\displaystyle
AJ^{2}\int_{-\infty}^{\infty}\frac{d\nu}{2\pi}\int\frac{d^{3}q}{8\pi^{3}}\mathit{Im}G(\mathbf{k+q+Q},\omega+\nu)$
(13)
$\displaystyle\times\mathit{Im}\chi(\mathbf{q+Q,}\nu)[f(\nu+\omega)+b(\nu)]$

where $f(\omega),b(\nu)$ are the Fermi and Bose functions respectively, and
$J$ is the coupling constant of electrons and spin fluctuations. The
conduction electron spectral function is approximated by

$\mathit{Im}G(\mathbf{k},\omega)=\pi\delta(\omega-\epsilon_{\mathbf{k}})$ (14)

We will employ an isotropic model of the conduction band with energy
$\epsilon_{k}=k^{2}/2m-\epsilon_{F}$, for which the angular integral may be
done, assuming $q\ll Q$, with the result

$\frac{1}{2}\int_{-1}^{1}d\cos\theta\mathit{Im}G(\mathbf{k+q+Q},\omega+\nu)\approx\frac{\pi}{2}\frac{m}{k_{F}Q}\approx\frac{1}{\epsilon_{F}}$
(15)

At low frequency $\omega<T$ we may drop $\omega$ in the argument of the Fermi
function, whence $\mathit{Im}\Sigma(\mathbf{k,}\omega)$ is approximately
independent of frequency

$\displaystyle\mathit{Im}\Sigma(\mathbf{k,}\omega$ $\displaystyle\lesssim$
$\displaystyle T)\approx
A\frac{n_{S}J^{2}}{I_{z}\epsilon_{F}}\int_{-\infty}^{\infty}\frac{d\nu}{2\pi\sinh(\nu/T)}$
$\displaystyle\times\int_{0}^{\infty}\frac{dqq^{2}}{2\pi^{2}}\frac{\nu/\gamma}{(\tau+(q/k_{F})^{2})^{2}+(\nu/\gamma)^{2}}$

The divergence of the momentum integral at $q\rightarrow 0$ is cut off in two
ways: (i) in the limit $\tau\rightarrow 0$ the Landau damping term provides
the cutoff just as is the case at the quantum critical point (QCP). The
ensuing quantum critical behavior is therefore found to extend in the phase
diagram from the QCP along the phase boundary in a narrow strip of width
$\tau_{x}T_{c}$, where $\tau_{x}\approx T_{c}/\gamma$; (ii) in the limit
$\nu\rightarrow 0$, or more generally $T\rightarrow 0$ (considering that the
frequency integral is confined to $|\nu|\lesssim T$) the cutoff is provided by
the term $\tau$ . This is the quantum disordered regime, in which
$\mathit{Im}\Sigma$ is found to depend critically on $\tau$. The latter
behavior may be seen in analogy to the approach to the QCP at $T=0$ from the
quantum disordered side.

Provided the transition temperature is not too high, $T_{c}\ll\gamma$ , there
exists a wide regime of reduced temperatures $\tau_{x}\ll\tau\ll 1$ with
$\tau_{x}\approx T_{c}/\gamma$ for which $\mathit{Im}\Sigma$ is approximately
given by

$\displaystyle\mathit{Im}\Sigma(\mathbf{k,}\omega$ $\displaystyle\lesssim
T)\approx\frac{3}{2\pi}A\frac{n_{S}nJ^{2}T}{I_{z}\epsilon_{F}\gamma}\int_{0}^{T}d\nu\int_{\sqrt{\tau}}\frac{dq/k_{F}}{(q/k_{F})^{2}}$
$\displaystyle\approx\frac{3}{2\pi}A\frac{n_{S}nJ^{2}}{I_{z}}\frac{T^{2}}{\epsilon_{F}\gamma}\frac{1}{\sqrt{\tau}}$
(17)

where $n=k_{F}^{3}/3\pi^{2}$ is the density of conduction electrons. The
resistivity follows as

$\rho_{sfl}^{3d}=\rho_{0}\frac{1}{\sqrt{\tau}},\text{ \ \
}\tau_{x}\ll\tau\text{ \ }$ (18)

where

$\rho_{0}\approx\frac{9}{4}R_{Q}k_{F}^{-1}A\frac{n_{S}J}{I_{z}}\frac{nJ}{\epsilon_{F}}\frac{T_{c}^{2}}{\epsilon_{F}\gamma},\text{
\ }$ (19)

and $R_{Q}=\frac{h}{e^{2}}\approx 25.81$k$\Omega$ is the quantum resistance.
In the ordered phase, in the temperature regime $\tau_{x,l}\ll|\tau|\ll 1$ ,
where $\tau_{x,l}\approx T_{c}/\gamma_{l}$, with $\gamma_{l}=\gamma
c_{l}M_{0}^{2}$, the resistivity scale $\rho_{0}$ in Eq. (18) is replaced by
$\rho_{0,l}=\rho_{0}/(M_{0}\sqrt{c_{l}})$. In the temperature regime close to
the transition, defined by $-\tau_{x,l}<\tau<\tau_{x}$ , the resistivity is
given by

$\rho_{sfl}^{3d}=\rho_{0}\frac{1}{\sqrt{\tau_{x}}},\text{ \ \ \
}-\tau_{x,l}<\tau<\tau_{x}$ (20)

We now estimate the $\rho_{sfl}$ of IrMn layers as studied in Ziman20 . The
reference resistance setting the scale is given by $R_{Q}k_{F}^{-1}\approx
258\mu\Omega$cm taking $k_{F}\approx 10^{8}$cm-1. The remaining factors at the
transition temperature are $\gamma\approx 10^{4}$K, $T_{c}/\gamma\approx 0.03$
, $I_{z}\approx T_{c}\approx 300$K, $nJ\approx N_{F}J\epsilon_{F}\approx
1.5\times 10^{3}$K, assuming $\epsilon_{F}\approx 10^{4}$K, and $A\approx 1$ ,
resulting in $\rho_{0}\approx 0.4\mu\Omega$cm. This is very small in
comparison to the observed resistivity (Fig.2 of the Supplementary information
of Tu et al. Ziman20  where $\rho\approx 170\mu\Omega$cm at $T_{c}$ ). Indeed,
the resistivity data do not show any trace of a peak at $T_{c}$. Here we
anticipate that the crossover scale $\tau_{x}\approx 0.06$ (see below),
leading to the estimate of the maximum height of the spin fluctuation induced
contribution as $\rho_{sfl}(T_{c})=\rho_{0}\frac{1}{\sqrt{\tau_{x}}}\approx
1.5\mu\Omega$cm. The model calculations for bulk MnIr Umetsu04  show a hump in
the electrical resistivity below $T_{c}$ .

###  III.2 Two-dimensional spin fluctuations

In this case the self-energy at low frequency is found as

$\displaystyle\mathit{Im}\Sigma(\mathbf{k,}\omega$ $\displaystyle\lesssim
T)\approx\frac{A}{2\pi^{2}}\frac{n_{S}J^{2}}{I_{z}\gamma}\int_{-\infty}^{\infty}\frac{d\nu\nu/\gamma}{\sinh(\nu/T)}$
$\displaystyle\times\int\frac{dqq/t_{AFM}}{(\tau+(q/k_{F})^{2})^{2}+(\nu/\gamma)^{2}}$
(21)
$\displaystyle\approx\frac{3}{2}A\frac{J^{2}nn_{S}}{I_{z}}\frac{T^{2}}{\gamma\epsilon_{F}}\frac{1}{k_{F}t_{AFM}}\frac{1}{\tau},$
(22)

using $n=k_{F}^{3}/3\pi^{2}$. We note that magnetic order is not destroyed by
transverse spin fluctuations because those are gapped out. The resistivity
shows a more strongly divergent behavior for $\tau\rightarrow 0$, up to the
crossover temperature $\tau_{x}$

$\rho_{sfl}^{2d}=\rho_{0}\frac{\pi}{k_{F}t_{AFM}}\frac{1}{\tau}\left\\{\begin{array}[]{c}1,\text{
\ \ \ \ \ \ \ \ \ }\tau\gg\tau_{x}\\\ \frac{1}{c_{l}M_{0}^{2}},\text{\ \
}-\tau\gg\tau_{x,l}\end{array}\right.$ (23)

The maximum of $\rho_{sfl}^{2d}$ is reached for $\tau<\tau_{x}$ and may be
estimated as
$\rho_{sfl}^{2d}\approx\rho_{0}\frac{\pi}{k_{F}t_{AFM}}\frac{1}{\tau_{x}}\approx
0.5\mu\Omega$cm, assuming $t_{AFM}=30$Å, which is again much less than the
total resistivity.

The system showing a record thermopower studied in Ziman20  consists of layers
of IrMn of thickness varying between $0.6-4$nm and is therefore quasi two-
dimensional. The momentum component $q_{z}$ of the spin fluctuations normal to
the sample plane is quantized, $q_{z,n}=\frac{\pi}{t_{AFM}}(2n+1)$,
$n=0,1,2...,$ where $t_{AFM}$ is the sample thickness. For sufficiently small
thickness $t_{AFM}$ only the lowest transverse mode is occupied. For this to
be valid we should have $(\frac{2\pi}{t_{AFM}}\xi_{0})^{2}\gg\tau_{x}$, where
$\xi_{0}\approx k_{F}^{-1}$ is the microscopic spin interaction length. For
the IrMn layers of thickness $\approx 3$nm the latter condition is not
satisfied, taking $\xi_{0}\approx 0.17$nm. We therefore conclude that the
three-dimensional model is more appropriate for describing the sample with
$t_{AFM}=2.8$nm, for which a detailed comparison with the thermopower data is
possible (see below).

##  IV Thermopower

The thermopower is defined by

$S=-\frac{\pi^{2}k_{B}^{2}T}{e}\frac{\rho^{\prime}(\epsilon_{F})}{\rho(\epsilon_{F})}$
(24)

where $\rho^{\prime}(E_{F})=d\rho/d\epsilon_{F}\approx$
$d\rho_{sfl}/d\epsilon_{F}$, assuming that of all contributions to the
resistivity, by scattering off phonons, off impurities, off magnetic ions, the
largest contribution is coming from critical spin fluctuations. The essential
dependence of $\rho_{sfl}$ on $\epsilon_{F}$ is through the shift of the
transition temperature $T_{c}$ induced by changing $\epsilon_{F}$ . Using
$d\tau/d\epsilon_{F}=-(T/T_{c}^{2})dT_{c}/d\epsilon_{F}$. one then finds above
the transition in case of three- or two-dimensional spin fluctuations

$S=\frac{\pi^{2}k_{B}^{2}}{e}\frac{\rho_{0}}{\rho}\frac{dT_{c}}{d\epsilon_{F}}(\frac{T}{T_{c}})^{2}\Theta(\tau-
T_{c}/\gamma)\left\\{\begin{array}[]{c}\tau^{-3/2}\text{, \ \ \ \ }d=3\\\
\frac{\pi}{k_{F}t_{AFM}}\tau^{-2}\text{,\ \ }d=2\end{array}\right.,$ (25)

a power law divergence $S\propto(T-T_{c})^{-\alpha}$, which is cut off at
$\tau\approx T_{c}/\epsilon_{F}$. The dependence of $T_{c}$ on $\epsilon_{F}$
has been discussed above,
$dT_{c}/d\epsilon_{F}\approx(n_{S}J/I_{z})(JN_{F})(T_{c}/\epsilon_{F})\approx(N_{F}J)^{2}\approx
0.02$. The peak height follows, using these estimates and $k_{B}/e\approx
8.6\times 10^{-3}$V/K , as $S(\tau_{x})\approx 240\mu$V/K, taking
$\tau_{x}=0.06$, for both three- or two-dimensional fluctuations, which is of
the order of magnitude observed in experiment, at least for the compounds
containing only IrMn magnetic layers. The compounds consisting of additional
adjacent CoFeB magnetic layers show even stronger thermopower, possibly
because the coupling constants $J,I$ are effectively changed by the magnetic
environment (increased $J$, decreased $I_{z}$). For temperatures below the
transition the prefactor $\rho_{0}$ and the cutoff scale $\tau_{x}$ are
replaced by $\rho_{0,l}$ and $\tau_{x,l}$.

Figure 2: Seebeck voltage of a IrMn nanostructure of thickness $t_{AFM}=2.8$nm
versus temperature. Theory according to Eq. (26) (solid line). Data points
from Ziman20 . The transition temperature to the antiferromagnetic phase is
fixed at $T_{c}=281$K.

In Fig. 2 we show a comparison of our theory with the Seebeck voltage data of
a sample with IrMn layer thickness of $t_{AFM}=2.8$nm. The transition
temperature is chosen to be $T_{c}=281$K, slightly less than the $285$K
estimated in Ziman20 . We use an interpolation expression connecting the
quantum critical regime at $|\tau|<\tau_{x}$ and the quantum disordered regime
at $|\tau|>\tau_{x}$ at $T>T_{c}$:
$|\tau|\Theta(\tau_{x}-|\tau|)+\tau_{x}\Theta(\tau_{x}-|\tau|)\approx\sqrt{\tau^{2}+\tau_{x}^{2}}$
(here $\Theta(x)$ is the unit step function). A similar expression is employed
for $T<T_{c}$. We find that the best fit is obtained assuming three-
dimensional spin fluctuations ($d=3$). The Seebeck voltage $V_{S}$ at
$T\gtrless T_{c}$ is then described by

$V_{S}^{>,<}=c_{1}^{>,<}(\frac{T}{T_{c}})^{2}[(\frac{T}{T_{c}}-1)^{2}+(\frac{T_{x}^{>,<}}{T_{c}}-1)^{2}]^{-0.75}+c_{2}^{>,<}$
(26)

The parameters used for fitting Eq. (26) to the data of Fig.2 of Ziman20  are
$T_{x}^{>}=299$K, $T_{x}^{<}=269.5$K, $c_{1}^{>}=0.89$mV, $c_{1}^{<}=0.4$mV,
$c_{2}^{>}=-6.5$mV, and $c_{2}^{<}=0$mV. The constants $c_{2}^{>,<}$ account
for contributions derived from other scattering processes (phonons,
impurities, magnetic ions) or from $d\rho_{0}/d\epsilon_{F}$.

##  V Conclusion

Usually the spin fluctuations in a spin system near a continuous thermal phase
transition into a magnetically ordered state are classical. In a metal,
however, the dynamics may have quantum character even at finite temperature on
account of coupling of the localized spins to the conduction electron spins,
giving rise to Landau damping. Provided the transition temperature is not too
high, $T_{c}\ll\gamma$, where $\gamma$ is the characteristic energy scale of
the Landau damping, quantum fluctuations may give rise to quantum critical
behavior. As we have shown here, the electrical resistivity acquires a quantum
critical contribution following a power law divergence in the reduced
temperature $|\tau|^{-\alpha}$ up to a crossover scale $\tau_{x}\ll 1$. While
the corresponding peak structure in the resistivity may be small in comparison
to the dominant terms due to scattering by phonons, magnetic moments or
impurities, it gives rise to a prominent peak in the thermopower. This is
because the critical temperature $T_{c}$ is found to depend sufficiently
strongly on the Fermi energy. We compare our theoretical results with recent
experimental observations Ziman20  on a heterostructure containing the
antiferromagnetic metal IrMn in the form of layers of several nanometers
thickness, and find excellent agreement. The strong dependence of the critical
temperature on the layer thickness found in Ziman20  may also be explained by
our theory.

While we have focussed on the experiments on IrMn, the results suggest that
large values of thermopower may result in other metallic films from the
coupling to magnetic fluctuations. The current experimental results argue for
fluctuations that are still in the three-dimensional regime, both from the
form of Figure 1, where the reduction by finite size of $T_{c,I}^{\ast}$ from
the limit $T_{c,I}$ is modest compared to the cutoff $\tau_{x}\approx
T_{c}/\gamma$ for the films of interest, and from the fit to the thermopower
data in Fig. 2. Note that we have argued that the much greater reduction in
critical temperature from bulk IrMn is due to a change in the microscopic
parameters (e.g. $I_{z}$) for the layers rather than a finite-size effect, in
the language of critical phenomena. One might expect even stronger enhancement
at room temperature by effective two-dimensionality, provided films of some
material are thin enough that the relative reduction due to finite size, i.e.
as in Equation 6, is greater than the cut-off $\tau_{x}$. For a fixed critical
temperature smaller values of $\tau_{x}$ would require a larger Landau damping
parameter $\gamma$, which seems unlikely. The effect of quantum fluctuations
at a finite temperature antiferromagnetic transition has also been seen in
specific heat data for two heavy fermion metals, CeCu6-xAux HvL94  and
YbRh2Si2 Krellner09  as will be shown in upcoming work WS21 .

##  VI Acknowledgments

We are grateful to Sadamichi Maekawa who suggested the importance of magnetic
fluctuations to thermopower. PW acknowledges support by a Distinguished Senior
Fellowship of Karlsruhe Institute of Technology. TZ thanks the Reimei program
of the ASRC, JAEA, Tokai, Japan for support.

## References

  * (1) S. Tu, T. Ziman, G. Yu,C. Wan, J. Hu, H. Wu, H. Wang, M. Liu, C. Guo, J. Zhang, M. A. Cabero Z., Y. Zhang, P. Gao, S. Liu, D. Yu, X. Han, I. Hallsteinsen, D. A. Gilbert, P. Wölfle, K. L. Wang, J-P- Ansermet, S. Maekawa, and H. Yu, Nature Commun.11, 2023 (2020).
  * (2) P. G. De Gennes, and J. Friedel, J. Phys. Chem. Solids 4, 71 (1958).
  * (3) M. E. Fisher, and J. S. Langer, Phys. Rev. Lett. 20, 665 (1968).
  * (4) O. Entin-Wolman, G. Deutscher, and R. Orbach, Phys. Rev. B 14, 4015 (1976).
  * (5) H. Hasegawa, and T. Moriya, J. Phys. Soc. Jpn. 36, 1542 (1974).
  * (6) T. Moriya, Spin Fluctuations in Itinerant Electron Magnetism, (Springer, New York, 1985).
  * (7) K. Ueda, J. Phys. Soc. Jpn. 43, 1497 (1977).
  * (8) R. Y. Umetsu, M. Miyakawa, K. Fukamichi, and A. Sakuma, Phys. Rev. B 69, 104411 (2004).
  * (9) L. Frangou, S. Oyarzun, S. Auffret, L. Vila, S. Gambarelli, and V. Baltz, Phys. Rev. Lett. 116, 077203 (2016).
  * (10) R. Zhang and R. F. Willis, Phys. Rev. Lett. 86, 2665 (2001).
  * (11) T. Moriya, Spin Fluctuations in Itinerant Electron Magnetism, Springer Series in Solid-State Sciences (Springer Berlin Heidelberg, 2012).
  * (12) J. R. Schrieffer, X. G. Wen, and S. C. Zhang, Phys. Rev. B 39, 11663 (1989).
  * (13) A. V. Chubukov and D. M. Frenkel, Phys. Rev. B 46, 11884 (1992).
  * (14) S. Adler, Phys. Rev. 137 B, 1022, (1965).
  * (15) H. v. Löhneysen, T. Pietrus, G. Portisch, H. G. Schlager, A. Schröder, M. Sieck, and T. Trappmann, Phys. Rev. Lett. 72, 3262 (1994).
  * (16) C. Krellner, S. Hartmann, A. Pikul, N. Oeschler, J. G. Donath, C. Geibel, F. Steglich, and J. Wosnitza, Phys. Rev. Lett. 102, 196402 (2009).
  * (17) P. Wölfle, and J. Schmalian, unpublished.


KANAZAWA-22-01

Semi-annihilating dark matter coupled with Majorons

Takumi<EMAIL_ADDRESS>and Takashi
<EMAIL_ADDRESS>
1Institute for Theoretical Physics, Kanazawa University, Kanazawa 920-1192,
Japan
2Institute of Liberal Arts and Science, Kanazawa University, Kanazawa
920-1192, Japan

##  1 Introduction

Weakly interacting massive particles (WIMPs) are one of the prominent thermal
dark matter candidates and are being searched for by various experiments and
observations such as indirect detection, direct detection, and collider
production. However, no clear signal of WIMPs has been found so far, and the
resultant bounds on the interactions between WIMPs and normal matter are
getting stronger. In particular, recent direct detection experiments impose a
stringent limit on the elastic scattering cross section with nuclei/nucleons.
The XENON1T and PandaX-4T Collaborations set a limit of $4.1\times
10^{-47}~{}\mathrm{cm}^{2}$ at a dark matter mass of $30~{}\mathrm{GeV}$ [1]
and $3.8\times 10^{-47}~{}\mathrm{cm}^{2}$ at $40~{}\mathrm{GeV}$ [2],
respectively.

Even in such a situation, thermally produced WIMPs are still attractive and
well-motivated because the dark matter relic abundance can be determined
independently from what occurred in the early universe. Furthermore, the
correctness of the thermal history of the universe has been partially proven
by successful Big Bang Nucleosynthesis (BBN). One of the ideas by which these
strong limits are naturally avoided is to consider a velocity-dependent cross
section such as pseudo-Nambu–Goldstone dark matter [3] or pseudo-scalar
interacting fermionic dark matter [4]. This type of dark matter may be
detected through experiments and observations if the dark matter is boosted by
some mechanism, because the interaction itself is not very small.

Another method to avoid the strong constraint of direct detection experiments
is to make dark matter one step closer to the dark sector, namely semi-
annihilating dark matter where semi-annihilations include (anti-)dark matter
particles in the final state [5, 6]. The relic abundance of semi-annihilating
dark matter is also thermally produced via the freeze-out mechanism, similar
to WIMPs, whereas it exhibits some features different from the canonical
WIMPs. For example, semi-annihilating dark matter may improve small-scale
problems, such as core-cusp, too-big-to-fail, and missing satellite problems
[7], via the core formation of dark matter halos [8, 9]. The self-heating of
dark matter by semi-annihilations also affects the calculation of the thermal
relic abundance [10], and the evolution of cosmological perturbations [11]. In
addition, semi-annihilations produce (semi-)relativistic boosted dark matter
in galaxies and exhibit distinctive signatures not induced by canonical WIMPs
[12, 13].

On the other hand, the mechanism of small neutrino mass generation is still
unknown, although neutrino oscillation experiments have confirmed that
neutrinos are massive [14]. The canonical seesaw mechanism is the simplest
possibility for small neutrino mass generation [15, 16, 17]; however, the mass
scale of the right-handed neutrinos is undetermined. It may be correlated with
a spontaneous breaking of global or gauge symmetries.

In this study, we construct a model of semi-annihilating dark matter and
simultaneously generate neutrino masses based on a global $U(1)_{B-L}$
symmetry. The Nambu–Goldstone boson associated with the global $U(1)_{B-L}$
symmetry, the so-called Majoron, can naturally be light enough compared with
the other particles because of its nature. Thus, the semi-annihilation
$\chi\chi\to\overline{\chi}J$ can easily be dominant over the other
annihilation channels, where $\chi$ is the dark matter and $J$ is the Majoron.
We calculate the thermal relic abundance of dark matter and find the parameter
space that can improve the small-scale problems via the self-heating effect by
semi-annihilation taking into account the constraints of the BBN and
perturbative unitarity of the couplings. Furthermore, we show that neutrinos
are produced by the subsequent decay of the Majoron produced by the semi-
annihilation $\chi\chi\to\overline{\chi}J\to\overline{\chi}\nu\nu$, and these
neutrinos indicate a characteristic box-shaped spectrum that can be detected
in future neutrino large detectors such as Hyper-Kamiokande (HK) [18, 19].
Finally, for comparison, we also investigate the case in which the Majoron
mass is larger than the dark matter mass. In this case, the three-to-two dark
matter annihilations, such as $\chi\chi\chi\to\chi\overline{\chi}$, determine
the relic abundance of dark matter instead of the semi-annihilation
$\chi\chi\to\overline{\chi}J$.

##  2 The model

| $Q$ | $u^{c}$ | $d^{c}$ | $L$ | $e^{c}$ | $H$ | $\nu_{R}^{c}$ | $\Phi$ | $\chi$
---|---|---|---|---|---|---|---|---|---
$SU(3)_{c}$ | $\bm{3}$ | $\overline{\bm{3}}$ | $\overline{\bm{3}}$ | $\bm{1}$ | $\bm{1}$ | $\bm{1}$ | $\bm{1}$ | $\bm{1}$ | $\bm{1}$
$SU(2)_{L}$ | $\bm{2}$ | $\bm{1}$ | $\bm{1}$ | $\bm{2}$ | $\bm{1}$ | $\bm{2}$ | $\bm{1}$ | $\bm{1}$ | $\bm{1}$
$U(1)_{Y}$ | $1/6$ | $2/3$ | $-1/3$ | $-1/2$ | $1$ | $1/2$ | $0$ | $0$ | $0$
$U(1)_{B-L}$ | $1/3$ | $-1/3$ | $-1/3$ | $-1$ | $1$ | $0$ | $1$ | $2$ | $2/3$
Table 1: Particle contents and quantum charges of the model.

We propose a model of semi-annihilating dark matter that is dominantly coupled
with neutrinos via a Majoron. The Standard Model (SM) is extended with a
global $U(1)_{B-L}$ symmetry, and the particle contents and charge assignments
are listed in Table 1. The new particles other than the SM fields are three
right-handed neutrinos $\nu_{Ri}^{c}~{}(i=1,2,3)$, the singlet complex scalar
$\Phi=(\phi+iJ)/\sqrt{2}$ including a Majoron field $J$, and another singlet
complex scalar $\chi$ which will be identified as dark matter. The Lagrangian
of the new particles invariant under the symmetry is given by

$\displaystyle\mathcal{L}=-\left(y_{\nu}\tilde{H}\overline{L}\nu_{R}+\frac{y_{\Phi}}{2}\Phi\overline{\nu_{R}^{c}}\nu_{R}+\mathrm{h.c.}\right)-\mathcal{V}(H,\Phi,\chi),$
(2.1)

where the generation indices are omitted. The scalar potential
$\mathcal{V}(H,\Phi,\chi)$ is given by

$\displaystyle\mathcal{V}(H,\Phi,\chi)=$
$\displaystyle-\frac{\mu_{H}^{2}}{2}|H|^{2}-\frac{\mu_{\Phi}^{2}}{2}|\Phi|^{2}+\frac{\mu_{\chi}^{2}}{2}|\chi|^{2}+\frac{\lambda_{H}}{2}|H|^{4}+\frac{\lambda_{\Phi}}{2}|\Phi|^{4}+\frac{\lambda_{\chi}}{2}|\chi|^{4}$
$\displaystyle+\lambda_{H\Phi}|H|^{2}|\Phi|^{2}+\lambda_{H\chi}|H|^{2}|\chi|^{2}+\lambda_{\Phi\chi}|\Phi|^{2}|\chi|^{2}$
$\displaystyle+\left(\frac{\lambda}{3\sqrt{2}}\Phi^{*}\chi^{3}-\frac{m_{J}^{2}}{4}\Phi^{2}+\mathrm{h.c.}\right).$
(2.2)

The last term $m_{J}^{2}\Phi^{2}/4$ is the explicit global $U(1)_{B-L}$
breaking term, which is essential for giving a mass to the Nambu–Goldstone
boson (Majoron), as we show in the following. The CP phases in the couplings
$\lambda$ and $m_{J}^{2}$ can be absorbed by the field redefinitions of $\chi$
and $\Phi$, respectively. Thus, all couplings in the scalar potential can
generally be regarded as real parameters.

###  2.1 The scalar sector

The electroweak gauge symmetry and global $U(1)_{B-L}$ symmetry are
spontaneously broken by the vacuum expectation values (VEVs) of $H$ and
$\Phi$, respectively. Then the fields $H$ and $\Phi$ can be parametrized as

$\displaystyle H=\frac{1}{\sqrt{2}}\left(\begin{array}[]{c}0\\\
v+h\end{array}\right),\quad\Phi=\frac{1}{\sqrt{2}}\left(v_{\phi}+\phi+iJ\right),$
(2.5)

where $v$ and $v_{\phi}$ are the VEVs. The CP-even components $h$ and $\phi$
mix with each other, whereas the CP-odd component $J$ itself is the mass
eigenstate being the Nambu–Goldstone boson associated with the global
$U(1)_{B-L}$ symmetry, the so-called Majoron. Using the stationary conditions
given by

$\displaystyle\mu_{H}^{2}$
$\displaystyle=\lambda_{H}v^{2}+\lambda_{H\Phi}v_{\phi}^{2},$ (2.6)
$\displaystyle\mu_{\Phi}^{2}$
$\displaystyle=\lambda_{\Phi}v_{\phi}^{2}+\lambda_{H\Phi}v^{2}-m_{J}^{2},$
(2.7)

the mass matrix of the CP-even components is written as

$\displaystyle\mathcal{V}\supset\frac{1}{2}\left(\begin{array}[]{cc}h&\phi\end{array}\right)\left(\begin{array}[]{cc}\lambda_{H}v^{2}&\lambda_{H\Phi}vv_{\phi}\\\
\lambda_{H\Phi}vv_{\phi}&\lambda_{\Phi}v_{\phi}^{2}\end{array}\right)\left(\begin{array}[]{c}h\\\
\phi\end{array}\right).$ (2.13)

This mass matrix can be diagonalized by a unitary matrix. As a result, the
gauge eigenstates $h$ and $\phi$ are rewritten by the unitary matrix and mass
eigenstates $h_{1}$ and $h_{2}$ as

$\displaystyle\left(\begin{array}[]{c}h\\\
\phi\end{array}\right)=\left(\begin{array}[]{cc}\cos\theta&\sin\theta\\\
-\sin\theta&\cos\theta\end{array}\right)\left(\begin{array}[]{c}h_{1}\\\
h_{2}\end{array}\right),$ (2.20)

where $h_{1}$ is the SM-like Higgs boson with a mass of
$m_{h_{1}}=125~{}\mathrm{GeV}$, and $h_{2}$ is the additional Higgs boson with
mass $m_{h_{2}}$.

The Higgs mixing angle $\sin\theta$ is constrained by K and B meson decays,
BBN and supernova explosions, which are summarized in Ref. [20]. The upper
bound is roughly given by $\sin\theta\lesssim 3\times 10^{-4}$ for
$m_{h_{2}}\lesssim 5~{}\mathrm{GeV}$ via K and B meson decays. For
$m_{h_{2}}\lesssim 0.2~{}\mathrm{GeV}$, the lower bound $\sin\theta\gtrsim
10^{-5}$ is also imposed by BBN and supernova explosions. The Higgs invisible
decays $h_{1}\to\chi\chi,h_{2}h_{2}$ also give an upper bound on the Higgs
mixing angle. However this bound is not strong as that from K and B meson
decays.

The complex scalar $\chi$ can be stabilized by the remnant $\mathbb{Z}_{3}$
symmetry after spontaneous symmetry breaking. Thus, it can be a dark matter
candidate whose mass is given by

$\displaystyle
m_{\chi}^{2}=\frac{1}{2}\left(\mu_{\chi}^{2}+\lambda_{H\chi}v^{2}+\lambda_{\Phi\chi}v_{\phi}^{2}\right).$
(2.21)

###  2.2 The neutrino sector

The full neutrino mass matrix arises from the Yukawa couplings after
spontaneous symmetry breaking:

$\displaystyle\mathcal{L}_{\nu}^{\mathrm{mass}}=$
$\displaystyle-\overline{\nu_{L}}m_{D}\nu_{R}-\frac{1}{2}\overline{\nu_{R}^{c}}M_{R}\nu_{R}+\mathrm{h.c.}$
$\displaystyle=$
$\displaystyle~{}-\left(\begin{array}[]{cc}\overline{\nu_{L}}&\overline{\nu_{R}^{c}}\\\
\end{array}\right)\left(\begin{array}[]{cc}0&m_{D}\\\
m_{D}^{T}&M_{R}\end{array}\right)\left(\begin{array}[]{c}\nu_{L}^{c}\\\
\nu_{R}\end{array}\right)+\mathrm{h.c.}$ (2.27) $\displaystyle\equiv$
$\displaystyle~{}-\left(\begin{array}[]{cc}\overline{\nu_{L}}&\overline{\nu_{R}^{c}}\\\
\end{array}\right)M_{\nu}\left(\begin{array}[]{c}\nu_{L}^{c}\\\
\nu_{R}\end{array}\right)+\mathrm{h.c.},$ (2.31)

where $m_{D}\equiv y_{\nu}v/\sqrt{2}$ and $M_{R}\equiv
y_{\Phi}v_{\phi}/\sqrt{2}$. The above $6\times 6$ neutrino mass matrix
$M_{\nu}$ can be diagonalized by the unitary matrix $V$ as [21]

$\displaystyle V^{T}M_{\nu}V=\mathrm{diag}(m_{1},\cdots,m_{6}),$ (2.32)

and the gauge eigenstates of the neutrinos can be rewritten in terms of the
mass eigenstates $\nu_{i}$ as

$\displaystyle\left(\begin{array}[]{c}\nu_{L}^{c}\\\
\nu_{R}\end{array}\right)=V\left(\begin{array}[]{c}\nu_{1}^{c}\\\ \vdots\\\
\nu_{6}^{c}\end{array}\right),\hskip
28.45274pt\left(\begin{array}[]{c}\nu_{L}\\\
\nu_{R}^{c}\end{array}\right)=V^{*}\left(\begin{array}[]{c}\nu_{1}\\\
\vdots\\\ \nu_{6}\end{array}\right),$ (2.43)

where $\nu_{i}~{}(i=1,2,3)$ corresponds to the active neutrinos with small
neutrino masses of $\mathcal{O}(0.1)~{}\mathrm{eV}$, and $\nu_{i}~{}(i=4,5,6)$
are sterile neutrinos with heavy masses.

The interactions between the Majoron and neutrinos arise from the Yukawa
coupling $y_{\Phi}$. Using the neutrino mass eigenbasis, the interactions can
be written as [22, 23]

$\displaystyle\mathcal{L}_{\nu}^{\mathrm{int}}=$
$\displaystyle-\frac{iJ}{2v_{\phi}}\sum_{i=1}^{6}\sum_{j=1}^{6}\overline{\nu_{i}}\left[C_{ij}\left(m_{i}P_{L}-m_{j}P_{R}\right)+C_{ji}\left(m_{j}P_{L}-m_{i}P_{R}\right)+\frac{m_{i}+m_{j}}{2}\delta_{ij}\gamma_{5}\right]\nu_{j}$
$\displaystyle\approx$
$\displaystyle~{}\frac{iJ}{2v_{\phi}}\sum_{i=1}^{6}m_{i}\overline{\nu_{i}}\gamma_{5}\nu_{i},$
(2.44)

where the coefficient $C_{ij}$ is defined by
$C_{ij}\equiv\sum_{k=1}^{3}V_{ki}V_{kj}^{*}$, and satisfies the relation
$C_{ij}=C_{ji}^{*}$. Note that we have taken only the leading term
$C_{ij}\approx\delta_{ij}$ in the second line of Eq. (2.44). The Majoron $J$
decays into a pair of neutrinos $J\to\nu_{i}\nu_{i}~{}(i=1,2,3)$ via the
interactions, and the decay width $\Gamma_{J}$ is calculated as

$\displaystyle\Gamma_{J}=\sum_{i=1}^{3}\frac{m_{J}}{16\pi}\frac{m_{i}^{2}}{v_{\phi}^{2}}.$
(2.45)

The sum of the active neutrino masses is constrained by the cosmological bound
$\sum_{i=1}^{3}m_{i}<0.11~{}\mathrm{eV}$ [26]. In addition, the stronger limit

$\displaystyle\displaystyle 2.6\times
10^{-3}~{}\mathrm{eV}^{2}<\sum_{i=1}^{3}m_{i}^{2}<5.2\times
10^{-3}~{}\mathrm{eV}^{2}$ (2.46)

is imposed by the global fit of the neutrino oscillation data [14]. We take
the averaged value $3.9\times 10^{-3}~{}\mathrm{eV^{2}}$ in the following
calculations. Successful BBN may be affected if the lifetime of the Majoron is
longer than $1~{}\mathrm{s}$. This sets a bound on the parameters as

$\displaystyle\left(\frac{m_{J}}{100~{}\mathrm{MeV}}\right)\left(\frac{10~{}\mathrm{GeV}}{v_{\phi}}\right)^{2}\gtrsim
8.5.$ (2.47)

Additional couplings with quarks and charged leptons are induced at the one-
loop level [21, 22], and couplings with gauge bosons are induced at the two-
loop level [23]. However, these interactions are anticipated to be suppressed
sufficiently because of the loop factor and the small Dirac neutrino mass
squared $m_{D}m_{D}^{\dagger}$ in the current model. These couplings can be
important only if the sterile neutrino masses are heavy enough, such as the
canonical seesaw scale $M_{R}\sim 10^{14}~{}\mathrm{GeV}$ [22, 23].

##  3 Dark matter

###  3.1 Thermal relic abundance

The complex scalar $\chi$ can be dark matter stabilized by the remnant
$\mathbb{Z}_{3}$ symmetry. The interactions of the dark matter $\chi$ are
induced by the scalar potential. Considering the plausible mass spectrum
$m_{J}\ll m_{\chi}\ll m_{h_{1},h_{2}}$, the main dark matter annihilation
channels are $\chi\overline{\chi}\to JJ$ and the semi-annihilation
$\chi\chi\to\overline{\chi}J~{}(\overline{\chi}\overline{\chi}\to\chi J)$. The
annihilation cross sections for these channels are calculated as

$\displaystyle\sigma_{\chi\overline{\chi}}{v}_{\mathrm{rel}}=$
$\displaystyle~{}\frac{1}{16\pi
s}\sqrt{1-\frac{4m_{J}^{2}}{s}}\left[-\frac{\lambda_{H\chi}v}{2v_{\phi}}\sin\theta\cos\theta\left(\frac{m_{h_{1}}^{2}}{s-m_{h_{1}}^{2}}-\frac{m_{h_{2}}^{2}}{s-m_{h_{2}}^{2}}\right)\right.$
$\displaystyle~{}\hskip
88.2037pt-\frac{\lambda_{H\chi}v}{2v_{\phi}}\sin\theta\cos\theta\left(m_{h_{1}}^{2}-m_{h_{2}}^{2}\right)\left(\frac{\cos^{2}\theta}{s-m_{h_{1}}^{2}}+\frac{\sin^{2}\theta}{s-m_{h_{2}}^{2}}\right)$
$\displaystyle~{}\hskip
88.2037pt+\left.\lambda_{\Phi\chi}\left(1+\frac{m_{h_{1}}^{2}\sin^{2}\theta}{s-m_{h_{1}}^{2}}+\frac{m_{h_{2}}^{2}\cos^{2}\theta}{s-m_{h_{2}}^{2}}\right)\right]^{2},$
(3.1) $\displaystyle\sigma_{\chi\chi}{v}_{\mathrm{rel}}=$
$\displaystyle~{}\sigma_{\overline{\chi}\overline{\chi}}{v}_{\mathrm{rel}}=\frac{\lambda^{2}}{8\pi
s}\sqrt{1-\frac{\left(m_{\chi}+m_{J}\right)^{2}}{s}}\sqrt{1-\frac{\left(m_{\chi}-m_{J}\right)^{2}}{s}},$
(3.2)

where $s$ is the Mandelstam variable, and $v_{\mathrm{rel}}$ is the relative
velocity of the dark matter. In the above calculation, we included the factor
$1/2!$ for identical particles in the final state, but not in the initial
state in accordance with the standard calculations of quantum field theory
[24]. The relative sizes of these cross sections were controlled by the
quartic couplings $\lambda_{H\chi}$, $\lambda_{\Phi\chi}$, and $\lambda$. In
this study, we focus on the case of
$\lambda_{H\chi},\lambda_{\Phi\chi}\ll\lambda$ because we are interested in
the effects of semi-annihilations.

The number densities of dark matter and anti-dark matter follow the Boltzmann
equations:

$\displaystyle\frac{dn_{\chi}}{dt}+3Hn_{\chi}=$
$\displaystyle-\langle\sigma_{\chi\chi}{v}_{\mathrm{rel}}\rangle\left(n_{\chi}^{2}-n_{\overline{\chi}}n_{\chi}^{\mathrm{eq}}\right)+\frac{\langle\sigma_{\overline{\chi}\overline{\chi}}{v}_{\mathrm{rel}}\rangle}{2}\left(n_{\overline{\chi}}^{2}-n_{\chi}n_{\chi}^{\mathrm{eq}}\right),$
(3.3) $\displaystyle\frac{dn_{\overline{\chi}}}{dt}+3Hn_{\overline{\chi}}=$
$\displaystyle+\frac{\langle\sigma_{\chi\chi}{v}_{\mathrm{rel}}\rangle}{2}\left(n_{\chi}^{2}-n_{\overline{\chi}}n_{\chi}^{\mathrm{eq}}\right)-\langle\sigma_{\overline{\chi}\overline{\chi}}{v}_{\mathrm{rel}}\rangle\left(n_{\overline{\chi}}^{2}-n_{\overline{\chi}}n_{\chi}^{\mathrm{eq}}\right).$
(3.4)

Note that the Majoron $J$ is in thermal equilibrium throughout the evolution
of (anti-)dark matter number densities in the Boltzmann equations because of
the moderate Higgs mixing angle. Assuming no asymmetry between the dark matter
and the anti-dark matter particles ($n_{\chi}=n_{\overline{\chi}}\equiv
n/2$),111This assumption is reasonable because there is no CP violation in the
scalar potential. the total number density of the dark matter and anti-dark
matter particles defined by $n$ follows the combined Boltzmann equation:

$\displaystyle\frac{dn}{dt}+3Hn=$
$\displaystyle-\frac{\langle\sigma_{\chi\chi}{v}_{\mathrm{rel}}\rangle}{4}\left(n^{2}-nn^{\mathrm{eq}}\right).$
(3.5)

Note that the factor $1/4$ appears on the right-hand side compared with the
case of the canonical WIMPs because we consider the semi-annihilation of the
complex scalar dark matter here [25]. Although the Boltzmann equation is
numerically solved to reproduce the observed relic abundance $\Omega
h^{2}\approx 0.12$ in the later part of the paper [26], the semi-analytic
solution can also easily be obtained [6]. The semi-annihilation cross section
in Eq. (3.2) can be simplified as $\sigma_{\chi\chi}v_{\mathrm{rel}}\approx
3\lambda^{2}/(128\pi m_{\chi}^{2})$ in $m_{J}\ll m_{\chi}$ and the non-
relativistic limit. Then, because the typical magnitude of the required cross
section is $\sigma_{\chi\chi}v_{\mathrm{rel}}\sim 2\times
10^{-25}~{}\mathrm{cm^{3}/s}$ for the semi-annihilating complex scalar dark
matter, one can obtain the required magnitude of the quartic coupling
$\lambda$ as

$\displaystyle\lambda\approx 1.6\times
10^{-3}\left(\frac{m_{\chi}}{\mathrm{GeV}}\right),$ (3.6)

to reproduce the observed relic abundance.

###  3.2 Halo core formation via semi-annihilation

It has been proposed that semi-annihilating dark matter can form cores of dark
matter halos because the dark matter mass is converted into kinetic energy via
semi-annihilation, and the dark matter particles in the inner region of the
halos are thermalized with pressure [8]. This effect alleviates small-scale
problems such as core-cusp, too-big-to-fail, and missing satellite problems
[7]. Core formation via semi-annihilating dark matter occurs more frequently
for dwarf-sized halos than for more massive halos in the same time scale. This
is a different feature from core formation via strongly self-interacting dark
matter [7]. Another feature of the core formation via semi-annihilating dark
matter is dark matter velocity dependence, as we show in the following.

The order of the semi-annihilation cross section required for core formation
was estimated as [8]

$\displaystyle\langle\sigma_{\chi\chi}{v}_{\mathrm{rel}}\rangle=\frac{10^{-29}~{}\mathrm{cm^{3}/s}}{\xi}\left(\frac{m_{\chi}}{100~{}\mathrm{MeV}}\right)\left(\frac{\sigma_{0}}{5~{}\mathrm{km/s}}\right)^{2}\left(\frac{0.5~{}M_{\odot}~{}\mathrm{pc}^{-3}}{\rho_{c}}\right)\left(\frac{10^{10}~{}\mathrm{yr}}{t_{\mathrm{age}}}\right),$
(3.7)

where $\xi$ is the energy absorption efficiency given by the ratio between the
dark matter halo radius and the mean free path of the dark matter

$\displaystyle\xi=10^{-3}\left(\frac{r_{s}}{5~{}\mathrm{kpc}}\right)\left(\frac{\rho_{c}}{M_{\odot}~{}\mathrm{pc}^{-3}}\right)\left(\frac{\sigma_{\mathrm{self}}/m_{\chi}}{10^{-3}~{}\mathrm{cm^{2}/g}}\right),$
(3.8)

$\sigma_{0}$ is the dark matter velocity dispersion due to thermalization,
$\rho_{c}$ is the dark matter core density of the halo, $r_{s}$ is the halo
radius, $M_{\odot}$ is the solar mass, $t_{\mathrm{age}}$ is the age of the
universe, and $\sigma_{\mathrm{self}}$ is the self-interacting cross section
of dark matter. Note that the order of the self-interacting cross section
required for core formation is typically $\sigma_{\mathrm{self}}/m_{\chi}\sim
10^{-3}~{}\mathrm{cm^{2}/g}$, which is much smaller than the case of strongly
self-interacting dark matter requiring $\mathcal{O}(1)~{}\mathrm{cm^{2}/g}$
[7].

In the current model, the self-interacting cross sections can be calculated as
[27]

$\displaystyle\sigma_{\chi\chi\to\chi\chi}=$
$\displaystyle~{}\sigma_{\overline{\chi}\overline{\chi}\to\overline{\chi}\overline{\chi}}=\frac{1}{32\pi
m_{\chi}^{2}}\left(\lambda_{\chi}+\frac{\lambda^{2}v_{\phi}^{2}}{6m_{\chi}^{2}}\right)^{2},$
(3.9) $\displaystyle\sigma_{\chi\overline{\chi}\to\chi\overline{\chi}}=$
$\displaystyle~{}\frac{1}{16\pi
m_{\chi}^{2}}\left(\lambda_{\chi}-\frac{\lambda^{2}v_{\phi}^{2}}{2m_{\chi}^{2}}\right)^{2},$
(3.10)

in non-relativistic limit. Although additional contributions are induced by
the couplings $\lambda_{H\chi}$ and $\lambda_{\Phi\chi}$, these contributions
are subdominant owing to the current setup
$\lambda_{H\chi},\lambda_{\Phi\chi}\ll\lambda$. Because there is no asymmetry
between the dark matter and the anti-dark matter particles, the effective
self-interacting cross section can be defined by

$\displaystyle\sigma_{\mathrm{self}}^{\mathrm{eff}}\equiv$
$\displaystyle~{}\frac{1}{4}\sigma_{\chi\chi\to\chi\chi}+\frac{1}{4}\sigma_{\chi\overline{\chi}\to\chi\overline{\chi}}+\frac{1}{4}\sigma_{\overline{\chi}\overline{\chi}\to\overline{\chi}\overline{\chi}}$
$\displaystyle=$ $\displaystyle~{}\frac{1}{64\pi
m_{\chi}^{2}}\left[\left(\lambda_{\chi}+\frac{\lambda^{2}v_{\phi}^{2}}{6m_{\chi}^{2}}\right)^{2}+\left(\lambda_{\chi}-\frac{\lambda^{2}v_{\phi}^{2}}{2m_{\chi}^{2}}\right)^{2}\right],$
(3.11)

and quantity $\sigma_{\mathrm{self}}/m_{\chi}$ in Eq. (3.8) should be replaced
by the effective cross section
$\sigma_{\mathrm{self}}^{\mathrm{eff}}/m_{\chi}$.

###  3.3 Perturbative unitarity bound

The perturbativity bounds for the quartic couplings $\lambda$ and
$\lambda_{\chi}$ are simply given by $\lambda,\lambda_{\chi}<4\pi$. The
unitarity bounds can also be imposed by the calculation of the scattering
amplitudes for $\chi\chi\to\chi\chi$ and
$\chi\overline{\chi}\to\chi\overline{\chi}$, which are given by [28, 29]

$\displaystyle\sqrt{1-\frac{4m_{\chi}^{2}}{s}}\left[2\lambda_{\chi}+\frac{\lambda^{2}v_{\phi}^{2}}{s-m_{\chi}^{2}}\right]$
$\displaystyle\leq 16\pi,$ (3.12)
$\displaystyle\sqrt{1-\frac{4m_{\chi}^{2}}{s}}\left[2\lambda_{\chi}-\frac{\lambda^{2}v_{\phi}^{2}}{s-4m_{\chi}^{2}}\log\left(\frac{s-3m_{\chi}^{2}}{m_{\chi}^{2}}\right)\right]$
$\displaystyle\leq 8\pi,$ (3.13)

respectively. In the high-energy limit $s\to\infty$, the bounds are simplified
to the usual ones $\lambda_{\chi}\leq 8\pi$ and $\lambda_{\chi}\leq 4\pi$ and,
thus, no unitarity bound for the coupling $\lambda$. However, the above
unitarity bounds in Eqs. (3.12) and (3.13) should be satisfied in all ranges
of the Mandelstam variable $s$. Thus, we vary $s$ in the range of
$4m_{\chi}^{2}\lesssim s<\infty$ to obtain a conservative bound for the
coupling $\lambda$. The results are presented in Fig. 1, and we find that the
bound $\lambda v_{\phi}/m_{\chi}\lesssim 17.5$ is numerically obtained from
the figure.

Figure 1: Unitarity bounds obtained from the scattering processes
$\chi\chi\to\chi\chi$ and $\chi\overline{\chi}\to\chi\overline{\chi}$. Colored
regions are excluded.

###  3.4 Direct detection

The dark matter $\chi$ can scatter a nucleon off via the couplings
$\lambda_{H\chi}$ and $\lambda_{\Phi\chi}$ in the scalar potential. The
constraint of dark matter direct detection is not so severe because the dark
matter mass scale is sub-GeV in our setup as we will see later. However for
completeness, we evaluate the upper bound on these couplings below. The
elastic scattering cross section between dark matter and nucleon at zero
momentum transfer is calculated as

$\displaystyle\sigma_{\chi
N}=\frac{m_{N}^{4}v_{\phi}^{4}f_{N}^{2}\left(\lambda_{H\chi}\lambda_{\Phi}-\lambda_{\Phi\chi}\lambda_{H\Phi}\right)^{2}}{4\pi(m_{\chi}+m_{N})^{2}m_{h_{1}}^{4}m_{h_{2}}^{4}},$
(3.14)

where $f_{N}=\sum_{q}f_{N}^{q}$ is the scalar quark form factors, which can be
evaluated as $f_{p}=$0.284 for a proton and $f_{n}=0.287$ for a nucleon [30].
Taking the upper bound on the elastic cross section
$\sigma_{\mathrm{el}}^{\mathrm{exp}}\sim 10^{-40}~{}\mathrm{cm}^{2}$ at
$m_{\chi}=1~{}\mathrm{GeV}$ [31], we can derive the bound

$\displaystyle\frac{v_{\phi}^{4}}{m_{h_{2}}^{4}}\left(\lambda_{H\chi}\lambda_{\Phi}-\lambda_{\Phi\chi}\lambda_{H\Phi}\right)^{2}\lesssim
0.011.$ (3.15)

one can see that this bound is not strong because of the relation
$v_{\phi}\lesssim m_{h_{2}}$. The future direct detection experiment SuperCDMS
is sensitive to sub-GeV scale dark matter and is anticipated to update the
bound up to $10^{-43}~{}\mathrm{cm}^{2}$ [32].

###  3.5 Box-shaped neutrino spectra

Figure 2: Diagram of the neutrino production from the semi-annihilation.

The semi-annihilation processes $\chi\chi\to\overline{\chi}J$
($\overline{\chi}\overline{\chi}\to\chi J$) produce the on-shell Majoron $J$.
Then, the neutrinos are produced by the subsequent decay of the Majoron
$J\to\nu_{i}\nu_{i}~{}(i=1,2,3)$, as shown in Fig. 2. These neutrinos can be a
signal of semi-annihilating dark matter in the model. The energy of the
(anti-)dark matter particle in the final state and the on-shell Majoron are
kinematically fixed by the energy–momentum conservation as

$\displaystyle
E_{\overline{\chi}}=\frac{5m_{\chi}^{2}-m_{J}^{2}}{4m_{\chi}},\quad
E_{J}=\frac{3m_{\chi}^{2}+m_{J}^{2}}{4m_{\chi}},$ (3.16)

because the dark matter particles in the initial state are non-relativistic.
The differential decay width of the neutrinos can be computed as [33]

$\displaystyle\frac{d\Gamma_{J}}{dE_{\nu}}=\frac{\Gamma_{J}}{\sqrt{E_{J}^{2}-m_{J}^{2}}}\theta\big{(}E_{\nu}-E_{\nu}^{\mathrm{min}}\big{)}\theta\big{(}E_{\nu}^{\mathrm{max}}-E_{\nu}\big{)},$
(3.17)

where the minimum and maximum neutrino energy $E_{\nu}^{\mathrm{min}}$ and
$E_{\nu}^{\mathrm{max}}$ are given by

$\displaystyle
E_{\nu}^{\mathrm{min}}=\frac{E_{J}}{2}\left(1-\sqrt{1-\frac{m_{J}^{2}}{E_{J}^{2}}}\right),\quad
E_{\nu}^{\mathrm{max}}=\frac{E_{J}}{2}\left(1+\sqrt{1-\frac{m_{J}^{2}}{E_{J}^{2}}}\right),$
(3.18)

and the total decay width of Majoron $\Gamma_{J}$ is given by Eq. (2.45).
Therefore, the produced neutrino energy spectrum becomes box-shaped, as shown
by the purple lines in Fig. 3, where the neutrino energy spectrum
$dN_{\nu}/dx$ is defined by

$\displaystyle\frac{dN_{\nu}}{dx}\equiv\frac{2}{\Gamma_{J}}\frac{d\Gamma_{J}}{dx},$
(3.19)

with the dimensionless parameter $x\equiv E_{\nu}/m_{\chi}$. The factor $2$
comes from the fact that two neutrinos are produced for each Majoron decay. In
Fig. 3, the energy resolution of the experimental detectors $\sigma$ was
considered by [33]

$\displaystyle\frac{dN_{\nu}}{dx}=\int\frac{dN_{\nu}}{dx^{\prime}}\frac{1}{\sqrt{2\pi}\sigma}e^{-(x^{\prime}-x)^{2}/2\sigma^{2}}dx^{\prime}.$
(3.20)

Figure 3: Energy spectra of the neutrino $\nu$ produced by the semi-
annihilation processes where the energy resolution of $\sigma/x=10\%$ is taken
into account. The purple lines represent the spectra for the on-shell Majoron
$J$ ($m_{J}\ll m_{\chi}$), whereas the green lines represent the spectra for
the off-shell Majoron $J$ ($m_{\chi}<m_{J}$). See Section 3.6 for the off-
shell case.

In the case of semi-annihilating dark matter, the neutrino flux produced by
dark matter semi-annihilations in the Milky Way halo is given by [34]

$\displaystyle\frac{d\Phi_{\nu}}{dx}=\frac{\langle\sigma_{\chi\chi}v_{\mathrm{rel}}\rangle}{4}\mathcal{J}_{\mathrm{avg}}\frac{R_{0}\rho_{0}^{2}}{m_{\chi}^{2}}\frac{dN_{\nu}}{dx},$
(3.21)

where $\mathcal{J}_{\mathrm{avg}}$ is the $J$-factor averaged over the Milky
Way halo, $R_{0}=8.5~{}\mathrm{kpc}$ is the distance between the galactic
center and the solar system, and $\rho_{0}=0.3~{}\mathrm{GeV/cm^{3}}$ is the
local dark matter density in the solar system. Although the $J$-factor
generally depends on the chosen dark matter halo profiles, it is not sensitive
because the $J$-factor is averaged over a large region of the solid angle
[35]. Using data from Super-Kamiokande (SK) and the future sensitivity of HK,
the current bound and future prospects on the standard dark matter
annihilation cross section for the channel $\chi\chi\to\nu\nu$ in the mass
range $10~{}\mathrm{MeV}\lesssim m_{\chi}\lesssim 200~{}\mathrm{MeV}$ have
been studied with the averaged $J$-factor $\mathcal{J}_{\mathrm{avg}}=5$ in
[34, 36]. These bound and future prospects for the monochromatic neutrino
lines can be translated into the case of the box-shaped spectrum in our case
by comparing the predicted box-shaped neutrino flux at
$E_{\nu}=E_{\nu}^{\mathrm{max}}$, given by Eq. (3.18) to the monochromatic
flux at $E_{\nu}=m_{\chi}$. This is plausible because the predicted neutrino
flux at $E_{\nu}=E_{\nu}^{\mathrm{max}}$ is expected to provide the strongest
limit on the cross section in our model.

###  3.6 Off-shell Majoron case

Here, we consider the case of $m_{\chi}<m_{J}\ll m_{h_{2}}$ for completeness,
although it may be less motivated from the viewpoint of the light pseudo-
Nambu–Goldstone boson. In this case, the semi-annihilation
$\chi\chi\to\overline{\chi}J$ and its CP conjugate processes are kinematically
forbidden. However, three-body semi-annihilation
$\chi\chi\to\overline{\chi}\nu\nu$ is possible via the off-shell Majoron. The
differential cross section is calculated as

$\displaystyle\frac{d\sigma_{\chi\chi}{v}_{\mathrm{rel}}}{dx}=\frac{\lambda^{2}\Gamma_{J}}{16\pi^{2}m_{\chi}^{3}m_{J}}\left[\frac{z}{1-z}+\log(1-z)\right],$
(3.22)

where the dimensionless parameter $z$ is defined as $z\equiv
wx\left(3-4x\right)/(1-x)$ with $w=m_{\chi}^{2}/m_{J}^{2}<1$. The total cross
section can be obtained by integrating Eq. (3.22) in the range $0<x<3/4$,

$\displaystyle\sigma_{\chi\chi}{v}_{\mathrm{rel}}=$
$\displaystyle~{}\frac{\lambda^{2}m_{J}\Gamma_{J}}{64\pi^{2}m_{\chi}^{4}}\Bigg{[}\frac{2-15w+9w^{2}}{2\sqrt{(1-w)(1-9w)}}\log\left(\frac{5-9w-3\sqrt{(1-w)(1-9w)}}{5-9w+3\sqrt{(1-w)(1-9w)}}\right)$
$\displaystyle~{}\hskip 48.36958pt-6w+(2-5w)\log{2}\Bigg{]}.$ (3.23)

Based on the above calculation, the neutrino spectrum defined by

$\displaystyle\frac{dN_{\nu}}{dx}\equiv\frac{2}{\sigma_{\chi\chi}v_{\mathrm{rel}}}\frac{d\sigma_{\chi\chi}v_{\mathrm{rel}}}{dx},$
(3.24)

is shown as the green lines in Fig. 3. The spectra can be significantly larger
than the box-shaped spectra at high energies. However, the bound and future
prospects for the cross section are much weaker than the on-shell Majoron case
studied above because the cross section in Eq. (3.23) is suppressed by the
small decay width of the Majoron $\Gamma_{J}$ given by Eq. (2.45).
Furthermore, this suppressed cross section was too small to reproduce the
observed relic abundance.222The other ways to detect the neutrinos for the
off-shell case at the DUNE experiment and colliders have been studied though
these cannot be simply applied to our case [37, 38].

Therefore, the relic abundance of dark matter is determined by three-to-two
self-annihilation processes such as $\chi\chi\chi\to\chi\overline{\chi}$ and
$\chi\overline{\chi}\chi\to\overline{\chi}\overline{\chi}$, instead of semi-
annihilations.333Another three-to-two annihilation channel
$\chi\overline{\chi}\chi\to\chi J$ opens if $m_{\chi}\lesssim m_{J}\lesssim
2m_{\chi}$. Assuming no CP asymmetry in the dark sector, as in the previous
case, the Boltzmann equation for the total dark matter number density $n$ is
given by [27]

$\displaystyle\frac{dn}{dt}+3Hn=-\frac{1}{4}\langle\overline{\sigma}_{3\to
2}v_{\mathrm{rel}}^{2}\rangle\left(n^{3}-n^{2}n_{\mathrm{eq}}\right),$ (3.25)

where $\langle\overline{\sigma}_{3\to 2}v_{\mathrm{rel}}^{2}\rangle$ is the
effective three-to-two cross section defined by

$\displaystyle\langle\overline{\sigma}_{3\to
2}v_{\mathrm{rel}}^{2}\rangle\equiv$
$\displaystyle~{}\frac{1}{3!}\langle\sigma_{\chi\chi\chi}v_{\mathrm{rel}}^{2}\rangle+\frac{1}{2!}\langle\sigma_{\chi\overline{\chi}\chi}v_{\mathrm{rel}}^{2}\rangle$
$\displaystyle=$
$\displaystyle~{}\frac{\sqrt{5}\lambda^{2}v_{\phi}^{2}}{1152\pi
m_{\chi}^{7}}\left[\frac{1}{3!}\left(\frac{\lambda^{2}v_{\phi}^{2}}{2m_{\chi}^{2}}+\lambda_{\chi}\right)^{2}+\frac{1}{2!}\frac{1}{48}\left(\frac{13\lambda^{2}v_{\phi}^{2}}{2m_{\chi}^{2}}-37\lambda_{\chi}\right)^{2}\right].$
(3.26)

In Eq. (3.26), the thermally averaged cross sections
$\langle\sigma_{\chi\chi\chi}v_{\mathrm{rel}}^{2}\rangle$ and
$\langle\sigma_{\chi\overline{\chi}\chi}v_{\mathrm{rel}}^{2}\rangle$ include
the symmetry factors for identical particles in the final state, but not in
the initial state as same as Eqs. (3.1) and (3.2), and only the $s$-wave
contribution is taken into account because the dark matter particles in the
initial state are non-relativistic.

The analytic solution for the observed relic abundance is approximately given
by $m_{\chi}\approx 0.035\alpha_{\mathrm{eff}}$, where the effective coupling
$\alpha_{\mathrm{eff}}$ is defined by $\langle\overline{\sigma}_{3\to
2}v_{\mathrm{rel}}^{2}\rangle\equiv\alpha_{\mathrm{eff}}^{3}/m_{\chi}^{5}$
[39, 40, 27]. Thus, for the simple case with $\lambda_{\chi}=0$, because the
effective coupling is evaluated as $\alpha_{\mathrm{eff}}\approx
0.034\lambda^{2}v_{\phi}^{2}/m_{\chi}^{2}$, we obtain the solution
$\lambda\approx
29\left(m_{\chi}/\mathrm{GeV}\right)^{1/2}\left(m_{\chi}/v_{\phi}\right)$.

##  4 Numerical results

We perform numerical computation and find the parameter space that can
simultaneously reproduce the observed dark matter relic abundance, form dark
matter halo cores and exhibit a neutrino signal from the Majoron decay. We
have five parameters relevant to the numerical computation: $m_{\chi}$,
$m_{J}$, $v_{\phi}$, $\lambda$, and $\lambda_{\chi}$. The Majoron mass should
be much smaller than the dark matter mass so that sufficient kinetic energy is
provided to the (anti-)dark matter particles by the semi-annihilation
processes, and the neutrino energy spectrum becomes box-shaped. Here we take
$m_{J}/m_{\chi}=0.1$ as a sample point. Note that the precise value of the
Majoron mass does not practically affect the relic abundance and neutrino
signals as long as $m_{J}/m_{\chi}\ll 1$, whereas the decay width of the
Majoron can be affected by Eq. (2.45).

The numerical result is shown in Fig. 4. In the orange region on the left-hand
side ($m_{\chi}\lesssim 20~{}\mathrm{MeV}$), the dark matter $\chi$ is still
in a thermal bath at the temperature $\sim 1~{}\mathrm{MeV}$ and affects the
effective number of neutrino species $N_{\mathrm{eff}}$ [41, 42]. In the green
region on the right-hand side, the successful BBN is spoiled by the additional
energy injection to the SM sector through $J\to\nu\nu$ because the Majoron is
too long-lived ($\tau_{J}>1~{}\mathrm{s}$), as shown in Eq. (2.45). The
relevant parameter space depends on the VEV $v_{\phi}$, and we have taken
$v_{\phi}/m_{\chi}=5,10$ in Fig. 4. The straight line whose color gradually
changes from blue to red corresponds to the parameter space that can reproduce
the observed relic abundance by semi-annihilation. The color represents the
value of the quartic coupling $\lambda_{\chi}$ required for the core formation
of dark matter halos. The perturbative unitarity bound for the quartic
coupling is imposed $\lambda_{\chi}<4\pi$. From Fig. 4, we can see that the
dark matter mass is $m_{\chi}\lesssim\mathcal{O}(1)~{}\mathrm{GeV}$ to
reproduce the observed relic abundance and realize the core formation of dark
matter halos via semi-annihilation simultaneously.

Figure 4: Parameter space in the ($m_{\chi},\lambda$) plane with
$m_{J}/m_{\chi}=0.1$. The gradient line from blue, black, and red represents
the parameter space that can reproduce the observed dark matter relic
abundance with the semi-annihilation $\chi\chi\to\overline{\chi}J$ where the
value of the quartic coupling $\lambda_{\chi}$ required for dark matter core
formation changes from $10^{-3}$ (blue) to $4\pi$ (red). The orange region on
the left-hand side is ruled out by $N_{\mathrm{eff}}$ (freeze-out temperature
$\lesssim 1~{}\mathrm{MeV}$). The green regions are excluded by the lifetime
of the Majoron ($\tau_{J}>1~{}\mathrm{s}$) for a fixed
$v_{\phi}/m_{\chi}=5,10$. The gray region is also excluded by the SK data
[34]. The dashed and dotted gray lines represent the HK sensitivity using Gd
10 years and without/with a boost factor of 10 (10BF) at the galaxies compared
with the value required for the relic abundance [36].

The gray region was excluded by the neutrino observation from the galactic
center by the SK [34]. The dashed gray line represents the future prospects of
the HK with gadolinium 10 years observation [36]. One can find that the
sensitivity of the HK can reach the parameter space favored for dark matter
relic abundance and core formation if the dark matter mass is in the range of
$20~{}\mathrm{MeV}\lesssim m_{\chi}\lesssim 35~{}\mathrm{MeV}$. The lower
dotted gray line is the future prospect with a boost factor of $10$ for the
semi-annihilation cross section. Such enhancement of the annihilation cross
section can be realized by mechanisms such as the Breit–Wigner effect [43, 44]
and the Sommerfeld effect [45, 46, 47, 48, 49]. For example, the simplest
extension of the model can be achieved by introducing a singlet scalar $S$
with a global $U(1)_{B-L}$ charge $-4/3$. Then, the semi-annihilation channel
$\chi\chi\to S^{*}\to\overline{\chi}J$ is enhanced at the resonance at
$2m_{\chi}\approx m_{S}$. If this kind of enhancement occurs, the dark matter
mass up to $m_{\chi}\sim 200~{}\mathrm{MeV}$ can be tested by HK with Gd, as
shown in Fig. 4.

Figure 5: Parameter space in the ($m_{\chi},\lambda$) plane where the Majoron
mass is fixed to be $m_{J}/m_{\chi}=3$. The purple region can reproduce the
observed relic abundance and form dark matter cores via strongly self-
interacting dark matter.
($0.1~{}\mathrm{cm^{2}/g}\leq\sigma_{\mathrm{self}}/m_{\chi}\leq
1~{}\mathrm{cm^{2}/g}$). The gray region is excluded by the unitarity bound
for the quartic coupling $\lambda$, and the other colored regions are the same
as in Fig. 4.

The numerical results for the off-shell Majoron case are shown in Fig. 5. The
purple region can reproduce the observed relic abundance and form dark matter
halo cores simultaneously. Here, we have taken the self-interacting cross
section required for the core formation as
$0.1~{}\mathrm{cm^{2}/g}\leq\sigma_{\mathrm{self}}/m_{\chi}\leq
1~{}\mathrm{cm^{2}/g}$ [7]. The upper gray region is excluded by the unitarity
bound for quartic coupling $\lambda$. The other colored regions were also
excluded for the same reasons as in Fig. 4. For $m_{\chi}\gtrsim
40~{}\mathrm{MeV}$, the purple region approximately coincides with the
analytical solution for the relic abundance $\lambda\approx
29\left(m_{\chi}/\mathrm{GeV}\right)^{1/2}\left(m_{\chi}/v_{\phi}\right)$. In
this mass region, the quartic coupling $\lambda$ is dominant over the other
coupling $\lambda_{\chi}$ to determine the relic abundance and form the dark
matter core formation. For the smaller mass region, because these couplings
are comparable or $\lambda_{\chi}$ is dominant, a wider region of $\lambda$
can be chosen at a fixed dark matter mass. Compared with the semi-annihilation
case ($m_{J}<m_{\chi}$), one finds that the quartic coupling $\lambda$ is
closer to the unitarity bound because a stronger self-interacting cross
section is required for dark matter halo core formation.

##  5 Summary and discussion

Although the canonical WIMPs are strongly constrained by recent experimental
and observational results, the thermal production mechanism of dark matter is
still attractive. Semi-annihilating dark matter is an alternative to thermal
dark matter, and semi-annihilations indicate some interesting phenomenological
aspects which do not emerge for the standard annihilations of dark matter,
such as core formation of dark matter halos, and boosted dark matter in the
final states. In this paper, we have proposed a model of semi-annihilating
dark matter based on global $U(1)_{B-L}$ symmetry, where two complex scalars
and three right-handed neutrinos are introduced as new fields to the SM. The
complex scalar is stabilized by the remnant $\mathbb{Z}_{3}$ symmetry after
spontaneous symmetry breaking. In addition, the small neutrino masses are also
generated via the seesaw mechanism. Because the pseudo-Nambu–Goldstone boson
associated with the global $U(1)_{B-L}$ symmetry, the so-called Majoron, is
anticipated to be light enough, the dark matter semi-annihilation into the
Majoron is naturally regarded as the dominant channel compared with the other
annihilation channels.

We have explored a parameter space that can reproduce the observed relic
abundance via the canonical freeze-out mechanism, form dark matter halo cores,
and be consistent with the cosmological observations and the perturbative
unitarity bounds. We have found that the dark matter mass should be in the
range of $20~{}\mathrm{MeV}\lesssim m_{\chi}\lesssim 1~{}\mathrm{GeV}$, and
the quartic coupling $\lambda$ that induces semi-annihilation should be in the
range of $3\times 10^{-5}\lesssim\lambda\lesssim 2\times 10^{-3}$.

In addition, the box-shaped spectrum of neutrinos can be generated by semi-
annihilation in the galaxies as a characteristic signature of the model. The
signal can be detected at the future large-volume neutrino detector, HK, if
the dark matter mass is rather light, $25~{}\mathrm{MeV}\lesssim
m_{\chi}\lesssim 35~{}\mathrm{MeV}$. Furthermore, a larger dark matter mass
region can be tested if the semi-annihilation cross section is enhanced by a
mechanism such as the Breit–Wigner and Sommerfeld effects.

For comparison, we have also investigated the case in which the Majoron mass
is heavier than the dark matter mass. In this case, the dark matter relic
abundance is determined by the three-to-two annihilation processes rather than
the semi-annihilation processes. As a result, stronger self-interactions are
required for relic abundance and dark matter core formation. In the off-shell
Majoron case, we have found that neutrinos are not generated as a
characteristic signal, unlike in the semi-annihilation case.

## Acknowledgments

This work was supported by a JSPS Grant-in-Aid for Scientific Research KAKENHI
Grant No. JP20K22349. The numerical computation in this work was carried out
at the Yukawa Institute Computer Facility.

## References

  * [1] E. Aprile et al. [XENON], Phys. Rev. Lett. 121, no.11, 111302 (2018) [arXiv:1805.12562 [astro-ph.CO]].
  * [2] Y. Meng et al. [PandaX-4T], [arXiv:2107.13438 [hep-ex]].
  * [3] C. Gross, O. Lebedev and T. Toma, Phys. Rev. Lett. 119, no.19, 191801 (2017) [arXiv:1708.02253 [hep-ph]].
  * [4] M. Freytsis and Z. Ligeti, Phys. Rev. D 83, 115009 (2011) [arXiv:1012.5317 [hep-ph]].
  * [5] T. Hambye, JHEP 01, 028 (2009) [arXiv:0811.0172 [hep-ph]].
  * [6] F. D’Eramo and J. Thaler, JHEP 06, 109 (2010) [arXiv:1003.5912 [hep-ph]].
  * [7] S. Tulin and H. B. Yu, Phys. Rept. 730, 1-57 (2018) [arXiv:1705.02358 [hep-ph]].
  * [8] X. Chu and C. Garcia-Cely, JCAP 07, 013 (2018) [arXiv:1803.09762 [hep-ph]].
  * [9] A. Kamada and H. J. Kim, Phys. Rev. D 102, no.4, 043009 (2020) [arXiv:1911.09717 [hep-ph]].
  * [10] A. Kamada, H. J. Kim, H. Kim and T. Sekiguchi, Phys. Rev. Lett. 120, no.13, 131802 (2018) [arXiv:1707.09238 [hep-ph]].
  * [11] A. Kamada, H. J. Kim and H. Kim, Phys. Rev. D 98, no.2, 023509 (2018) [arXiv:1805.05648 [hep-ph]].
  * [12] K. Agashe, Y. Cui, L. Necib and J. Thaler, JCAP 10, 062 (2014) [arXiv:1405.7370 [hep-ph]].
  * [13] T. Toma, [arXiv:2109.05911 [hep-ph]].
  * [14] P. F. de Salas, D. V. Forero, S. Gariazzo, P. Martínez-Miravé, O. Mena, C. A. Ternes, M. Tórtola and J. W. F. Valle, JHEP 02, 071 (2021) [arXiv:2006.11237 [hep-ph]].
  * [15] P. Minkowski, Phys. Lett. B 67, 421-428 (1977)
  * [16] T. Yanagida, Conf. Proc. C 7902131, 95-99 (1979) KEK-79-18-95.
  * [17] M. Gell-Mann, P. Ramond and R. Slansky, Conf. Proc. C 790927, 315-321 (1979) [arXiv:1306.4669 [hep-th]].
  * [18] K. Abe, T. Abe, H. Aihara, Y. Fukuda, Y. Hayato, K. Huang, A. K. Ichikawa, M. Ikeda, K. Inoue and H. Ishino, et al. [arXiv:1109.3262 [hep-ex]].
  * [19] N. F. Bell, M. J. Dolan and S. Robles, JCAP 09, 019 (2020) [arXiv:2005.01950 [hep-ph]].
  * [20] M. W. Winkler, Phys. Rev. D 99, no.1, 015018 (2019) [arXiv:1809.01876 [hep-ph]].
  * [21] A. Pilaftsis, Phys. Rev. D 49, 2398-2404 (1994) [arXiv:hep-ph/9308258 [hep-ph]].
  * [22] C. Garcia-Cely and J. Heeck, JHEP 05, 102 (2017) [arXiv:1701.07209 [hep-ph]].
  * [23] J. Heeck and H. H. Patel, Phys. Rev. D 100, no.9, 095015 (2019) [arXiv:1909.02029 [hep-ph]].
  * [24] M. E. Peskin and D. V. Schroeder, “An Introduction to quantum field theory”, Addison-Wesley (1995).
  * [25] P. Gondolo and G. Gelmini, Nucl. Phys. B 360, 145-179 (1991)
  * [26] N. Aghanim et al. [Planck], Astron. Astrophys. 641, A6 (2020) [erratum: Astron. Astrophys. 652, C4 (2021)] [arXiv:1807.06209 [astro-ph.CO]].
  * [27] S. M. Choi and H. M. Lee, JHEP 09, 063 (2015) [arXiv:1505.00960 [hep-ph]].
  * [28] L. Di Luzio, J. F. Kamenik and M. Nardecchia, Eur. Phys. J. C 77, no.1, 30 (2017) [arXiv:1604.05746 [hep-ph]].
  * [29] M. D. Goodsell and F. Staub, Eur. Phys. J. C 78, no.8, 649 (2018) [arXiv:1805.07306 [hep-ph]].
  * [30] G. Bélanger, F. Boudjema, A. Goudelis, A. Pukhov and B. Zaldivar, Comput. Phys. Commun. 231, 173-186 (2018) [arXiv:1801.03509 [hep-ph]].
  * [31] P. A. Zyla et al. [Particle Data Group], PTEP 2020, no.8, 083C01 (2020)
  * [32] R. Agnese et al. [SuperCDMS], Phys. Rev. D 95, no.8, 082002 (2017) [arXiv:1610.00006 [physics.ins-det]].
  * [33] A. Ibarra, S. Lopez Gehler and M. Pato, JCAP 07, 043 (2012) [arXiv:1205.0007 [hep-ph]].
  * [34] A. Olivares-Del Campo, C. Bœhm, S. Palomares-Ruiz and S. Pascoli, Phys. Rev. D 97, no.7, 075039 (2018) [arXiv:1711.05283 [hep-ph]].
  * [35] H. Yuksel, S. Horiuchi, J. F. Beacom and S. Ando, Phys. Rev. D 76, 123506 (2007) [arXiv:0707.0196 [astro-ph]].
  * [36] A. Olivares-Del Campo, S. Palomares-Ruiz and S. Pascoli, [arXiv:1805.09830 [hep-ph]].
  * [37] K. J. Kelly and Y. Zhang, Phys. Rev. D 99, no.5, 055034 (2019) [arXiv:1901.01259 [hep-ph]].
  * [38] K. J. Kelly, F. Kling, D. Tuckler and Y. Zhang, [arXiv:2111.05868 [hep-ph]].
  * [39] Y. Hochberg, E. Kuflik, T. Volansky and J. G. Wacker, Phys. Rev. Lett. 113, 171301 (2014) [arXiv:1402.5143 [hep-ph]].
  * [40] Y. Hochberg, E. Kuflik, H. Murayama, T. Volansky and J. G. Wacker, Phys. Rev. Lett. 115, no.2, 021301 (2015) [arXiv:1411.3727 [hep-ph]].
  * [41] C. Boehm, M. J. Dolan and C. McCabe, JCAP 12, 027 (2012) [arXiv:1207.0497 [astro-ph.CO]].
  * [42] C. Boehm, M. J. Dolan and C. McCabe, JCAP 08, 041 (2013) [arXiv:1303.6270 [hep-ph]].
  * [43] K. Griest and D. Seckel, Phys. Rev. D 43, 3191-3203 (1991)
  * [44] M. Ibe, H. Murayama and T. T. Yanagida, Phys. Rev. D 79, 095009 (2009) [arXiv:0812.0072 [hep-ph]].
  * [45] J. Hisano, S. Matsumoto and M. M. Nojiri, Phys. Rev. D 67, 075014 (2003) [arXiv:hep-ph/0212022 [hep-ph]].
  * [46] J. Hisano, S. Matsumoto and M. M. Nojiri, Phys. Rev. Lett. 92, 031303 (2004) [arXiv:hep-ph/0307216 [hep-ph]].
  * [47] J. Hisano, S. Matsumoto, M. M. Nojiri and O. Saito, Phys. Rev. D 71, 063528 (2005) [arXiv:hep-ph/0412403 [hep-ph]].
  * [48] J. Hisano, S. Matsumoto, O. Saito and M. Senami, Phys. Rev. D 73, 055004 (2006) [arXiv:hep-ph/0511118 [hep-ph]].
  * [49] J. Hisano, S. Matsumoto, M. Nagai, O. Saito and M. Senami, Phys. Lett. B 646, 34-38 (2007) [arXiv:hep-ph/0610249 [hep-ph]].