#  $D^{+}\to K_{s}^{0}\pi^{+}\eta$ reaction and $a_{0}(980)^{+}$

Natsumi Ikeno<EMAIL_ADDRESS>Department of Agricultural, Life and
Environmental Sciences, Tottori University, Tottori 680-8551, Japan  Jorgivan
M. Dias<EMAIL_ADDRESS>CAS Key Laboratory of Theoretical Physics,
Institute of Theoretical Physics, Chinese Academy of Sciences, Beijing 100190,
China  Wei-Hong Liang<EMAIL_ADDRESS>Department of Physics, Guangxi
Normal University, Guilin 541004, China  Guangxi Key Laboratory of Nuclear
Physics and Technology, Guangxi Normal University, Guilin 541004, China
Eulogio Oset<EMAIL_ADDRESS>Departamento de Física Teórica and IFIC,
Centro Mixto Universidad de Valencia-CSIC Institutos de Investigación de
Paterna, Aptdo.22085, 46071 Valencia, Spain  Department of Physics, Guangxi
Normal University, Guilin 541004, China

###### Abstract

We study the $D^{+}\to\bar{K}^{0}\pi^{+}\eta$ reaction where the $a_{0}(980)$
excitation plays a dominant role. We consider mechanisms of external and
internal emission at the quark level, hadronize the $q\bar{q}$ components into
two mesons and allow these mesons to undergo final state interaction where the
$a_{0}(980)$ state is generated. While the $a_{0}(980)$ production is the
dominant term, we also find other terms in the reaction that interfere with
this production mode and, through interference with it, lead to a shape of the
$a_{0}(980)$ significantly different from the one observed in other
experiments, with an apparently much larger width.

##  I Introduction

The $D^{+}\to K_{s}^{0}\pi^{+}\eta$ reaction measured by the BESIII
collaboration in Ref. BESIII:2020pxp , and more recently in Ref.
BESIII:2023htx  with more precision and an amplitude analysis, has turned into
an ideal reaction to isolate the $a_{0}(980)$ contribution. The reaction is
actually $D^{+}\to\bar{K}^{0}\pi^{+}\eta$, and the $\bar{K}^{0}$ is observed
as a $K_{s}^{0}$ state. It might look at a simple sight that this reaction is
just a copy of the $D^{0}\to K^{-}\pi^{+}\eta$ reaction measured by the Belle
collaboration Belle:2020fbd . Indeed, one goes from the first reaction to the
second by changing a $\bar{d}\to\bar{u}$ quark. Yet, the differences are
striking as we shall see in the present work, comparing it with the
theoretical study of the $D^{0}\to K^{-}\pi^{+}\eta$ reaction in Ref.
Toledo:2020zxj . The key point to see these striking differences is the fact
that the $K^{-}\pi^{+}$ in the $D^{0}\to K^{-}\pi^{+}\eta$ reaction can come
from $\bar{K}^{*0}$ excitation, but $\bar{K}^{0}\pi^{+}$ with positive charge
and an $s$ quark cannot come from $\bar{K}^{*0},K^{*-}$. Hence, there is no
$\bar{K}^{*}$ contribution in the $D^{+}\to\bar{K}^{0}\pi^{+}\eta$ reaction
which makes cleaner the $a_{0}(980)$ production as seen in the experiment
BESIII:2023htx .

The differences in the mass distributions in the $D^{0}\to K^{-}\pi^{+}\eta$
and $D^{+}\to\bar{K}^{0}\pi^{+}\eta$ reactions are striking. Indeed, in the
$D^{0}\to K^{-}\pi^{+}\eta$ reaction the $K^{-}\pi^{+}$ mass distribution is
dominated by the $\bar{K}^{*0}$ contribution, with a sharp peak at the
$\bar{K}^{*0}$ mass, while in the $D^{+}\to\bar{K}^{0}\pi^{+}\eta$ reaction
the $\bar{K}^{0}\pi^{+}$ mass distribution is rather structureless. This is
not all, because the $\pi^{+}\eta$ mass distribution in the $D^{0}\to
K^{-}\pi^{+}\eta$ reaction, shows indeed a peak for the $a_{0}(980)$
production but has large strength at low and high invariant masses which are
replicas of the $\bar{K}^{*0}$ resonance in the $K^{-}\pi^{+}$ channel. Even
more striking is the shape of the $K^{-}\eta$ mass distribution which has a
double hump shape created again by the presence of the $\bar{K}^{*0}$
resonance in the $K^{-}\pi^{+}$ channel. By contrast, the $\pi^{+}\eta$ mass
distribution in the $D^{+}\to\bar{K}^{0}\pi^{+}\eta$ reaction has a mass
distribution where the $a_{0}(980)$ dominates the spectrum and the
$\bar{K}^{0}\eta$ mass distribution has also very different shape.

It is worth mentioning that the related $D^{0}\to
K_{s}^{0}\pi^{+}\pi^{-},\bar{K}_{s}^{0}\pi^{0}\eta$ reactions were studied
prior to the experiments, paying attention to the $\pi^{+}\pi^{-}$ and
$\pi^{0}\eta$ mass distributions, predicting that a clear signal of the
$a_{0}(980)$ should be seen in these experiments Xie:2014tma .

Our purpose in the present work is to try to understand the spectrum from the
perspective that the $a_{0}(980)$ is a dynamically generated state from the
interaction of the $\pi\eta,K\bar{K}$ channels, which is well described by the
chiral unitary approach Oller:1997ti ; Kaiser:1998fi ; Markushin:2000fa ;
Nieves:1998hp . For this we shall investigate the reaction mechanism of
external and internal emission Chau:1982da  at the quark level, and see how
the hadronization of $q\bar{q}$ pairs can lead to the production of the final
state.

##  II Formalism

We study the $D^{+}\to\bar{K}^{0}\pi^{+}\eta$ reaction and look at the
mechanisms for external and internal emission.

###  II.1 External emission

In Fig. 1(a)

Figure 1: Mechanism with external emission at the quark level. (a)
$D^{+}\to\pi^{+}s\bar{d}$; (b) hadronization of the $s\bar{d}$ component
adding a $\bar{q}q$ pair producing two mesons.

we show the mechanism of external emission, Cabibbo favoured, and in Fig. 1(b)
the $s\bar{d}$ pair is hadronized in the following way

$s\bar{d}\to\sum_{i}s\;\bar{q}_{i}q_{i}\;\bar{d},$ (1)

which is easily interpreted writing the $q_{i}\bar{q}_{j}$ matrix in terms of
pseudoscalar mesons, $\mathcal{P}_{ij}$, which in the standard
$\eta-\eta^{\prime}$ mixing of Ref. Bramon:1992kr  is given by

$\mathcal{P}=\left(\begin{array}[]{cccc}\frac{1}{\sqrt{2}}\pi^{0}+\frac{1}{\sqrt{3}}\eta&\pi^{+}&K^{+}\\\\[5.69054pt]
\pi^{-}&-\frac{1}{\sqrt{2}}\pi^{0}+\frac{1}{\sqrt{3}}\eta&K^{0}\\\\[5.69054pt]
K^{-}&\bar{K}^{0}&~{}-\frac{1}{\sqrt{3}}\eta{}\\\ \end{array}\right).$ (2)

Hence,

$s\bar{d}\to\sum_{i}\;\mathcal{P}_{3i}\;\mathcal{P}_{i2}=\left(\mathcal{P}^{2}\right)_{32}=K^{-}\pi^{+}-\bar{K}^{0}\,\dfrac{\pi^{0}}{\sqrt{2}},$
(3)

where the $\bar{K}^{0}\eta$ component has cancelled.

The combination of Eq. (3), together with the $\pi^{+}$ of Fig. 1, does not
lead to the desired final state $\bar{K}^{0}\pi^{+}\eta$. However, through
rescattering the $K^{-}\pi^{+}$ and $\bar{K}^{0}\pi^{0}$ could lead to
$\bar{K}^{0}\eta$. Yet, the threshold of $\bar{K}^{0}\eta$ is $1040\textrm{
MeV}$, which is about $300\textrm{ MeV}$ above the peak of the
$K^{*}_{0}(700)$ (the kappa), where the $K^{-}\pi^{+}\to\bar{K}^{0}\eta$
amplitude has a reduced strength, and because of that we shall disregard this
contribution, as it was also done in Ref. Toledo:2020zxj .

Next, with the same mechanism of Fig. 1(a), we look at the hadronization of
the $\bar{d}u$ component, as shown in Fig. 2.

Figure 2: Hadronization of the $\bar{d}u$ component in external emission.

We have now

$\displaystyle u\bar{d}$ $\displaystyle\to$
$\displaystyle\sum_{i}u\;\bar{q}_{i}q_{i}\;\bar{d}=\sum_{i}\;\mathcal{P}_{1i}\;\mathcal{P}_{i2}=\left(\mathcal{P}^{2}\right)_{12}$
(4) $\displaystyle=$
$\displaystyle\dfrac{2}{\sqrt{3}}\,\eta\pi^{+}+K^{+}\bar{K}^{0},$

where now the $\pi^{0}\pi^{+}$ component has cancelled.

Hence we have now the combination

$H=\left(\frac{2}{\sqrt{3}}\,\eta\pi^{+}+K^{+}\bar{K}^{0}\right)\,\bar{K}^{0}.$
(5)

However, there is a subtlety here concerning the $K^{+}\bar{K}^{0}$ production
(see discussion in page 3 of Ref. Ikeno:2021kzf ) because for dynamical
reasons the $WK^{+}\bar{K}^{0}$ vertex goes as the difference of energies of
$K^{+}\bar{K}^{0}$ which vanishes in the average. However, due to the
different masses of $\eta\pi^{+}$ this cancellation does not occur, and
consequently we keep the $\eta\pi^{+}$ term and disregard the
$K^{+}\bar{K}^{0}$ one.

We can see that the term $H$ in Eq. (5) already contains the
$\bar{K}^{0}\pi^{+}\eta$ state at the tree level. In addition we can have
rescattering of the $\eta\pi^{+}$ and also of the $\bar{K}^{0}\pi^{+}$. Once
again, we disregard the $\eta\bar{K}^{0}$ rescattering for the reasons given
before. Then, the picture that we have for $\bar{K}^{0}\pi^{+}\eta$ production
is shown in Fig. 3, including tree level and rescattering.

Figure 3: Diagrams for $D^{+}\to\bar{K}^{0}\pi^{+}\eta$ coming from external
emission: (a) tree level; (b) $\eta\pi^{+}$ rescattering; (c)
$\bar{K}^{0}\pi^{+}$ rescattering.

Analytically we have for these terms originating from external emission

$\displaystyle t^{(ee)}$ $\displaystyle=$
$\displaystyle\mathcal{C}\left\\{h_{\eta\pi^{+}\bar{K}^{0}}\right.$
$\displaystyle+\left.h_{\eta\pi^{+}\bar{K}^{0}}\left[G_{\eta\pi^{+}}(M_{\rm
inv}(\eta\pi^{+}))\cdot t_{\eta\pi^{+},\,\eta\pi^{+}}(M_{\rm
inv}(\eta\pi^{+}))\right.\right.$
$\displaystyle+\left.\left.G_{\bar{K}^{0}\pi^{+}}(M_{\rm
inv}(\bar{K}^{0}\pi^{+}))\cdot
t_{\bar{K}^{0}\pi^{+},\,\bar{K}^{0}\pi^{+}}(M_{\rm
inv}(\bar{K}^{0}\pi^{+}))\right]\right\\},$

where $h_{\eta\pi^{+}\bar{K}^{0}}$ is the weight of the
$\eta\pi^{+}\bar{K}^{0}$ component in Eq. (5),

$h_{\eta\pi^{+}\bar{K}^{0}}=\dfrac{2}{\sqrt{3}},$ (7)

and $G_{i},t_{ij}$ are the loop functions of two mesons and the scattering
matrices for the transition of channel $i$ to channel $j$, taken as in Ref.
Toledo:2020zxj . $\mathcal{C}$ in Eq. (II.1) is a global constant that will be
used to get the normalization of the data.

###  II.2 Internal emission

For internal emission we have the diagrams of Fig. 4, which account for the
hadronization of quark pairs.

Figure 4: Diagrams for internal emission: (a) hadronization of the $s\bar{d}$
pair; (b) hadronization of the $u\bar{d}$ pair.

Following the path of subsection II.1, we have now for the hadronization

$\displaystyle s\bar{d}$ $\displaystyle\to$
$\displaystyle\sum_{i}s\;\bar{q}_{i}q_{i}\;\bar{d}=\sum_{i}\;\mathcal{P}_{3i}\;\mathcal{P}_{i2}=\left(\mathcal{P}^{2}\right)_{32}$
(8) $\displaystyle=$ $\displaystyle
K^{-}\pi^{+}-\dfrac{1}{\sqrt{2}}\,\pi^{0}\bar{K}^{0},$

where the $\bar{K}^{0}\eta$ terms have cancelled, and

$\displaystyle u\bar{d}$ $\displaystyle\to$
$\displaystyle\sum_{i}u\;\bar{q}_{i}q_{i}\;\bar{d}=\sum_{i}\;\mathcal{P}_{1i}\;\mathcal{P}_{i2}=\left(\mathcal{P}^{2}\right)_{12}$
(9) $\displaystyle=$
$\displaystyle\dfrac{2}{\sqrt{3}}\,\eta\pi^{+}+K^{+}\bar{K}^{0},$

where the $\pi^{+}\pi^{0}$ terms have cancelled. Summing the two terms we
have, including the $\pi^{+}$ in Fig. 4(a) and $\bar{K}^{0}$ in Fig. 4(b),

$H^{\prime}=K^{-}\pi^{+}\pi^{+}-\dfrac{1}{\sqrt{2}}\;\pi^{0}\pi^{+}\bar{K}^{0}+\dfrac{2}{\sqrt{3}}\;\eta\pi^{+}\bar{K}^{0}+K^{+}\bar{K}^{0}\bar{K}^{0}.$
(10)

Once again we have a tree level $\eta\pi^{+}\bar{K}^{0}$ production and the
other terms can lead to this final state through rescattering, as shown in
Fig. 5.

Figure 5: Diagrams for $\eta\pi^{+}\bar{K}^{0}$ production coming from
internal emission: (a) tree level; (b) $\pi^{+}\eta$ rescattering; (c)
$\pi^{+}\bar{K}^{0}$ rescattering; (d) $K^{+}\bar{K}^{0}$ rescattering.

In Fig. 5 we have disregarded the possible rescattering of
$K^{-}\pi^{+}\to\bar{K}^{0}\eta$, $\bar{K}^{0}\pi^{0}\to\bar{K}^{0}\eta$,
$\bar{K}^{0}\eta\to\bar{K}^{0}\eta$ for the same reasons as in the former
subsection. We have checked numerically that these terms are, indeed,
negligible.

Analytically we have for the diagrams of Fig. 5 corresponding to internal
emission

$\displaystyle t^{(ie)}$ $\displaystyle=$
$\displaystyle\beta\mathcal{C}\left\\{\bar{h}_{\eta\pi^{+}\bar{K}^{0}}\right.$
(11) $\displaystyle+$
$\displaystyle\left.\bar{h}_{\eta\pi^{+}\bar{K}^{0}}\left[G_{\eta\pi^{+}}(M_{\rm
inv}(\eta\pi^{+}))\cdot t_{\eta\pi^{+},\,\eta\pi^{+}}(M_{\rm
inv}(\eta\pi^{+}))\right.\right.$ $\displaystyle+$
$\displaystyle\left.G_{\pi^{+}\bar{K}^{0}}(M_{\rm
inv}(\pi^{+}\bar{K}^{0}))\cdot
t_{\pi^{+}\bar{K}^{0},\,\pi^{+}\bar{K}^{0}}(M_{\rm
inv}(\pi^{+}\bar{K}^{0}))\right]$ $\displaystyle+$ $\displaystyle
2\left.\bar{h}_{K^{+}\bar{K}^{0}\bar{K}^{0}}\;G_{K^{+}\bar{K}^{0}}(M_{\rm
inv}(\pi^{+}\eta))\right.$ $\displaystyle\left.\ \ \cdot
t_{K^{+}\bar{K}^{0},\pi^{+}\eta}(M_{\rm inv}(\pi^{+}\eta))\right\\},$

where the weights $\bar{h}_{i}$ are now

$\bar{h}_{\eta\pi^{+}\bar{K}^{0}}=\dfrac{2}{\sqrt{3}};~{}~{}~{}\bar{h}_{K^{+}\bar{K}^{0}\bar{K}^{0}}=1,$
(12)

and the factor $2$ multiplying $\bar{h}_{K^{+}\bar{K}^{0}\bar{K}^{0}}$ in Eq.
(11) accounts for the two $\bar{K}^{0}$ identical particles. In Eq. (11) we
have put a weight $\beta\mathcal{C}$, conserving the global factor
$\mathcal{C}$. The $\beta$ should be a factor of the order of
$\frac{1}{N_{c}}$ ($N_{c}$ number of colors) which one expects from the ratio
of terms from internal emission to those of external emission.

The total amplitude for the process is now given by

$t=t^{(ee)}+t^{(ie)}.$ (13)

###  II.3 $K^{*}_{0}(1430)$ contribution

In Ref. BESIII:2023htx  it was found that the scalar $K^{*}_{0}(1430)$
[$I\,(J^{P})=\frac{1}{2}\,(0^{+})$] state contributed to the process and shows
up in the $K\eta$ mass distribution. The process can proceed as shown in Fig.
6.

Figure 6: Diagrammatic representation of
$D^{+}\to\pi^{+}\bar{K}^{*}_{0}(1430)\to\pi^{+}\bar{K}^{0}\eta$.

We take into account this contribution phenomenologically by means of the
amplitude

$t^{*}=\mathcal{D}\,e^{i\phi}\;\dfrac{M_{D}^{2}}{s_{13}-M^{2}_{K^{*}_{0}}+iM_{K^{*}_{0}}\,\Gamma_{K^{*}_{0}}},$
(14)

where $\mathcal{D}$ and $\phi$ will be chosen as free parameters and
$s_{13}=(p_{\bar{K}^{0}}+p_{\eta})^{2}$ taking the order of the particles
$\bar{K}^{0}(1),\pi^{+}(2),\eta(3)$. The factor $M_{D}^{2}$ is put to have
$\mathcal{D}$ dimensionless. The final amplitude will now be

$t=t^{(ee)}+t^{(ie)}+t^{*}.$ (15)

Then in order to calculate the mass distributions, we use the PDG standard
formula pdg2023

$\dfrac{d^{2}\Gamma}{ds_{12}\;ds_{23}}=\dfrac{1}{(2\pi)^{3}}\;\dfrac{1}{32\,M^{3}_{D}}\;|t|^{2}.$
(16)

We can integrate over the limits of the PDG formula pdg2023  to get
$d\Gamma/ds_{12}$ integrating over $s_{23}$. Then, by cyclical permutations in
the formulas we can get $d\Gamma/ds_{13}$, $d\Gamma/ds_{23}$ and compare with
experiment BESIII:2023htx .

###  II.4 Scattering amplitudes

We need the amplitudes

$t_{\eta\pi^{+},\,\eta\pi^{+}},~{}~{}t_{\bar{K}^{0}\pi^{+},\,\bar{K}^{0}\pi^{+}},~{}~{}t_{\bar{K}^{0}K^{+},\,\pi^{+}\eta}.$

We use the results of Ref. Lin:2021isc  (Eq. (A.4) of Ref. Lin:2021isc ),
which considers explicitly the $\eta-\eta^{\prime}$ mixing, and for the
channels $K^{+}K^{-}(1),K^{0}\bar{K}^{0}(2),\pi^{0}\eta(3)$ one finds the
matrix elements of the potential

$\displaystyle V_{11}=-\dfrac{s}{2f^{2}},$ $\displaystyle
V_{12}=-\dfrac{s}{4f^{2}},$ $\displaystyle
V_{13}=-\dfrac{3s-2m^{2}_{K}-m^{2}_{\eta}}{3\sqrt{6}f^{2}},$ $\displaystyle
V_{22}=-\dfrac{s}{2f^{2}},$ $\displaystyle
V_{23}=-\dfrac{3s-2m^{2}_{K}-m^{2}_{\eta}}{3\sqrt{6}f^{2}},$ $\displaystyle
V_{33}=-\dfrac{2m^{2}_{\pi}}{3f^{2}}.$ (17)

By evaluating the coupled channel $T$ matrix

$T=[1-VG]^{-1}\,V,$ (18)

we find, considering that $|\pi^{+}\rangle=-|11\rangle$ of isospin,

$t_{\eta\pi^{+},\,\eta\pi^{+}}=t_{\eta\pi^{0},\,\eta\pi^{0}}.$ (19)

Then, taking into account that $K^{+}K^{-}$ in term of $|I,\,I_{3}\rangle$ is

$K^{+}K^{-}=-\left(\dfrac{1}{\sqrt{2}}|10\rangle+\dfrac{1}{\sqrt{2}}|00\rangle\right),$
(20)

we find that

$t_{\bar{K}^{0}K^{+},\,\pi^{+}\eta}=\sqrt{2}\;t_{K^{+}K^{-},\,\pi^{0}\eta}.$
(21)

We also need the $K\pi\to K\pi$ amplitudes, which we take from the Appendix of
Ref. Toledo:2020zxj , using the channels $\pi^{-}K^{+}(1),\pi^{0}K^{0}(2),\eta
K^{0}(3)$. Once again, using isospin coefficients and $C$ parity
$C\bar{K}^{0}\pi^{+}=K^{0}\pi^{-}$, we obtain

$t_{\bar{K}^{0}\pi^{+},\,\bar{K}^{0}\pi^{+}}=\dfrac{2}{3}\,T_{22}+\dfrac{1}{3}\,T_{11}+\dfrac{2\sqrt{2}}{3}\,T_{12}.$
(22)

Note that this amplitude has $I=3/2$ and hence does not contain the
$K^{*}_{0}(700)$.

##  III Results

We have four parameters, $\mathcal{C}$, $\beta$, $\mathcal{D}$, and
$\exp(i\phi)$, in our theoretical model. To determine these parameters, we
perform a best fit to the three mass distributions of $\bar{K}^{0}\pi^{+}$,
$\bar{K}^{0}\eta$, and $\pi^{+}\eta$ with the experimental data of Ref.
BESIII:2023htx . The parameter set obtained from the fit is shown in Table 1.
The parameter $\mathcal{C}$ provides a global normalization for external and
internal emission. The parameter $\beta$ gives the relative weight of the
internal emission mechanism to the external emission, while $\mathcal{D}$
quantifies the strength of the $K_{0}^{*}(1430)$ contribution. Additionally,
the phase $\exp(i\phi)$ corresponds to the interference between the
$K_{0}^{*}(1430)$ and other contributions, influencing the non-trivial shapes
of the various mass distributions under consideration. In Figs. 7–9, we
illustrate the mass distribution results obtained with the parameters from
Table 1. It is evident that our theoretical calculations closely replicate the
experimental data.

In Fig. 7, we can see a clear peak dominating the $\pi^{+}\eta$ spectrum
around $1.0$ GeV, corresponding to the $a_{0}(980)$ resonance, which, in our
model, is encoded in the $T$-matrix elements $t_{\eta\pi^{+},\,\eta\pi^{+}}$,
and $t_{\bar{K}^{0}K^{+},\,\eta\pi^{+}}$, defined in Eqs. (19) and (21),
respectively. In contrast, both the $\bar{K}^{0}\pi^{+}$ scattering term in
Eq. (22) and the contributions from $K_{0}^{*}(1430)$ have relatively small
strengths compared to the $a_{0}(980)$ resonance, particularly around $1.0$
GeV, where their corresponding strengths are 8-9 times smaller. Consequently,
we can claim that the peak observed in the experiment can be identified as the
$a_{0}(980)$ state. Furthermore, it is crucial to note that due to
interference with other contributions, particularly with the tree level in
Eqs. (II.1) and (11), the lineshape of the $a_{0}(980)$ is broader than those
observed in other reactions. Indeed, the apparent width of approximately $150$
MeV exceeds the range of 50-100 MeV as stated in the PDG pdg2023 , or the
width observed in the $\chi_{c1}\to\eta\pi^{+}\pi^{-}$ decay BESIII:2016tqo
(see also Ref. Liang:2016hmr ).

In Fig. 8, we reproduce a double hump structure in the $\bar{K}^{0}\eta$ mass
distribution. The $K^{*}_{0}(1430)$ spectrum has a peak structure at
$M_{\bar{K}^{0}\eta}=M_{K^{*}_{0}}=1425$ MeV with the width
$\Gamma_{K^{*}_{0}}=270$ MeV as seen in Eq. (14). The hump structure comes
from the interference between the $K^{*}_{0}(1430)$ and $a_{0}(980)$
contributions. Thus, the phase exp$(i\phi)$ in Table 1 is required to be
negative.

Fig. 9 shows no distinct peak structure in the $\bar{K}^{0}\pi^{+}$ mass
distribution, except for a broad bump around $1.25$ GeV. The
$\bar{K}^{0}\pi^{+}$ spectrum also exhibits a relatively large contribution,
with a discontinuity at $1.05$ GeV. This behavior arises from the cut-off mass
$M_{\textrm{cut}}$ applied to the contribution from the products $G\cdot t$ in
Eqs. (II.1) and (11), following the prescription discussed in Refs.
Debastiani:2016ayp ; Toledo:2020zxj . The introduction of this cut is
essential for extrapolating $Gt$ to high energies, given that the two-body
amplitudes from the chiral unitary approach are applicable up to about $1200$
MeV. Furthermore, it is worth noting that, as discussed in Refs.
Debastiani:2016ayp ; Toledo:2020zxj , the overall results exhibit minimal
dependence on the parameters associated with these cuts. In this study, we
adopt $M_{\textrm{cut}}=1050$ MeV and $\alpha=0.0037$ MeV-1. According to our
findings, the impact of $M_{\textrm{cut}}$ is noticeable only in the
$\bar{K}^{0}\pi^{+}$ distribution. Moreover, when exploring different values
for the $M_{\textrm{cut}}$, such as $1150$ MeV as shown in subsection III.1,
we observe no significant influence on the distribution lineshapes, except in
the $\bar{K}^{0}\pi^{+}$ spectrum, where the dip shifts from $1050$ MeV to
$1150$ MeV. Thus, we conclude that the dip in our model is not physical, and
we could have a smooth curve in that region. The relevant thing is that this
change in $M_{\rm cut}$ has little repercussion in the $\bar{K}^{0}\eta$ mass
distribution and in particular in the $\pi^{+}\eta$ distribution, where the
effects are negligible. The latter one is the most relevant channel in the
present work, where we are concerned about the effect of the $a_{0}(980)$ in
this reaction, and, hence, the conclusions that we draw about the effect of
the $a_{0}(980)$ are not affected by the mentioned uncertainties.

Table 1:  Values of the parameters from the fit.  Parameters |
---|---
$\mathcal{C}$ |  $486.90$
$\mathcal{D}$ |  $63.94$
$\beta$ |  $0.70$
$\phi$ |  $-2.16$ radians
Figure 7: The $\pi^{+}\eta$ mass distribution $d\Gamma/dM_{\pi^{+}\eta}$. The
total theoretical result is shown in the red line. The contributions of the
$a_{0}(980)$, the $\bar{K}^{0}\pi^{+}$ scattering terms, and the
$K^{*}_{0}(1430)$ are shown in the blue, green, and magenta lines,
respectively. The experimental data shown by points are taken from Ref.
BESIII:2023htx .  Figure 8: The $\bar{K}^{0}\eta$ mass distribution
$d\Gamma/dM_{\bar{K}^{0}\eta}$. The same as Fig. 7 but for the
$\bar{K}^{0}\eta$ mass distribution.  Figure 9: The $\bar{K}^{0}\pi^{+}$ mass
distribution $d\Gamma/dM_{\bar{K}^{0}\pi^{+}}$. The same as Fig. 7 but for the
$\bar{K}^{0}\pi^{+}$ mass distribution.

###  III.1 Theoretical uncertainties

In this subsection, we show the results of several tests to determine some
uncertainties of our theoretical models to the three mass distributions.

####  III.1.1 Effect of the cut mass $M_{\rm cut}$

Figure 10 shows the results for the considered distributions with
$M_{\textrm{cut}}$ fixed at $1150$ MeV. Following a new fit in this setting,
we obtain the parameters $\mathcal{C}=532.04$, $\mathcal{D}=90.28$,
$\beta=0.70$, and $\phi=-2.07$ radians. As a result, while there are
noticeable changes in the strengths of individual contributions, the
lineshapes for the first two distributions show only minimal differences,
especially in the $\pi^{+}\eta$ distribution where the $a_{0}(980)$ peak is
evident. Conversely, the dip in the $\bar{K}^{0}\pi^{+}$ spectrum has shifted
to $1150$ MeV compared to its position in Fig. 9. As discussed earlier, this
dip is directly influenced by the parameter $M_{\textrm{cut}}$. Nevertheless,
we observe only slight changes in the strength and shape of the
$\bar{K}^{0}\pi^{+}$ distribution.

Figure 10:  The mass distributions of $\pi^{+}\eta$ (left), $\bar{K}^{0}\eta$
(middle), and $\bar{K}^{0}\pi^{+}$ (right) with fixed $M_{\rm cut}=1150$ MeV.
The parameters $\mathcal{C}=532.04$, $\mathcal{D}=90.28$, $\beta=0.70$, and
$\phi=-2.07$ radians are used.

####  III.1.2 Effect of the parameter $\beta$

The parameter $\beta$ gives the relative weight of the internal emission
mechanism to the external emission. The value of the parameter $\beta$ might
be expected to be of the order of $1/N_{c}$. Here we restrict the value of the
$\beta$ within $[-0.33:0.33]$ and we make a fit. In Fig. 11, we show the three
mass distributions with $\beta=0.33$ which is obtained from the best fit. The
parameters obtained from the fit are $\mathcal{C}=691.80$,
$\mathcal{D}=71.29$, $\beta=0.33$, and $\phi=-2.29$ radians. Once again, we
see that the changes are not significant.

Figure 11:  The mass distributions of $\pi^{+}\eta$ (left), $\bar{K}^{0}\eta$
(middle), and $\bar{K}^{0}\pi^{+}$ (right) with fixed $M_{\rm cut}=1050$ MeV.
The parameters $\mathcal{C}=691.80$, $\mathcal{D}=71.29$, $\beta=0.33$, and
$\phi=-2.29$ radians are used.

####  III.1.3 Effect of the $K_{0}^{*}(1430)$ mass

The mass of the $K_{0}^{*}(1430)$ has a relatively large uncertainty,
$M_{K^{*}_{0}}=1425\pm 50$ MeV pdg2023 . Here we take $M_{K^{*}_{0}}=1385$
MeV, which is consistent within the error bars, and then we perform the same
calculations for the three invariant mass distributions. The parameters
obtained from the fit are $\mathcal{C}=473.34$, $\mathcal{D}=57.27$,
$\beta=0.70$, and $\phi=-2.39$ radians. We show the three mass distribution in
Fig. 12. The resulting calculation of the $\bar{K}^{0}\pi^{+}$ mass
distribution is in better agreement with the data because the peak position of
$K^{*}_{0}(1430)$ has moved a bit to the left from Fig. 8.

Figure 12:  The mass distributions of $\pi^{+}\eta$ (left), $\bar{K}^{0}\eta$
(middle), and $\bar{K}^{0}\pi^{+}$ (right) with fixed $M_{K^{*}_{0}}=1385$ MeV
and fixed $M_{\rm cut}=1050$ MeV. The parameters $\mathcal{C}=473.34$,
$\mathcal{D}=57.27$, $\beta=0.70$, and $\phi=-2.39$ radians are used.

In conclusion, we can see the clear peak of the $a_{0}(980)$ contribution even
considering the uncertainties. We found that the
$D^{+}\to\bar{K}^{0}\pi^{+}\eta$ reaction is a good reaction to see the peak
structure of $a_{0}(980)$ in the $\pi^{+}\eta$ mass distribution. Yet, we also
explained why the observed shape and width do not correspond to those seen in
other experiments, through the interference with other terms that appear
necessarily linked to the $a_{0}(980)$ production in our theoretical approach.

##  IV Conclusions

We have performed an analysis of the $D^{+}\to\bar{K}^{0}\pi^{+}\eta$ reaction
based on the picture of the $a_{0}(980)$ resonance as dynamically generated
from the interaction of the $\pi\eta$, $K\bar{K}$ channels, which emerges from
the study of the chiral unitary approach and has been tested in many previous
reactions. We show that this reaction is drastically different from the
apparently analogous one $D^{0}\to K^{-}\pi^{+}\eta$, and we trace it to the
absence of a $K^{*}$ contribution in the $D^{+}\to\bar{K}^{0}\pi^{+}\eta$
reaction, while it is the driving term in the $D^{0}\to K^{-}\pi^{+}\eta$ one.
This leads to a much cleaner signal of the $a_{0}(980)$ excitation in the
$D^{+}\to\bar{K}^{0}\pi^{+}\eta$ reaction as seen in the experiment.

In our study, we begin by looking at the dominant weak decay modes of external
and internal emission at the quark level. Then we proceed with the
hadronization of the $q\bar{q}$ pairs into two mesons and finally we allow the
meson pairs to interact. We also take into account the contribution of the
$\bar{K}^{*}_{0}(1430)$ as done in the experimental analysis and fit it to the
data. Up to this empirical information and a global normalization constant,
our framework only depends on the relative strength of the internal emission
mechanism to the external emission one, which is also a fit parameter, and we
get a result barely consistent with the expected large $N_{c}$ reduction. With
this framework, we obtain a fair reproduction of the three mass distributions
and we observe that the driving term in this reaction is the excitation of the
$a_{0}(980)$. Yet, we note that the shape obtained has a much larger width
than the one observed in other experiments, and we trace back this feature to
the interference of the dominant $a_{0}(980)$ excitation with the other
mechanisms that are also generated in the final state interaction of the
mesons produced, together with the $\bar{K}^{*}_{0}(1430)$ excitation. This
explanation is important to prevent taking the present shape of the
$a_{0}(980)$ as a measure of the width of the $a_{0}(980)$ state, and
reconcile the findings of the $D^{+}\to\bar{K}^{0}\pi^{+}\eta$ reaction with
the shapes of the $a_{0}(980)$ seen in other experiments.

###### Acknowledgements.

One of us, W.H.L would like to thank Prof. Bai-Cian Ke and Prof. Liao-Yuan
Dong for helpful discussion. N. I. and J. M. D. would like to express
gratitude to Guangxi Normal University for their warm hospitality, as part of
this work was conducted there. This work is partly supported by the National
Natural Science Foundation of China (NSFC) under Grants No. 11975083 and No.
12365019, and by the Central Government Guidance Funds for Local Scientific
and Technological Development, China (No. Guike ZY22096024). J. M. Dias
acknowledges the support from the Chinese Academy of Sciences under Grants No.
XDB34030000 and No. YSBR-101; by the National Key R&D Program of China under
Grant No. 2023YFA1606703. This project has received funding from the European
Union Horizon 2020 research and innovation programme under the program
H2020-INFRAIA-2018-1, grant agreement No. 824093 of the STRONG-2020 project.

## References

  * (1) M. Ablikim et al. [BESIII], Phys. Rev. Lett. 124, 241803 (2020).
  * (2) M. Ablikim et al. [BESIII], [arXiv:2309.05760 [hep-ex]].
  * (3) Y. Q. Chen et al. [Belle], Phys. Rev. D 102, 012002 (2020).
  * (4) G. Toledo, N. Ikeno and E. Oset, Eur. Phys. J. C 81, 268 (2021).
  * (5) J. J. Xie, L. R. Dai and E. Oset, Phys. Lett. B 742, 363-369 (2015).
  * (6) J. A. Oller and E. Oset, Nucl. Phys. A 620, 438-456 (1997); [erratum: Nucl. Phys. A 652, 407-409 (1999)]
  * (7) N. Kaiser, Eur. Phys. J. A 3, 307-309 (1998).
  * (8) V. E. Markushin, Eur. Phys. J. A 8, 389-399 (2000).
  * (9) J. Nieves and E. Ruiz Arriola, Phys. Lett. B 455, 30-38 (1999).
  * (10) L. L. Chau, Phys. Rept. 95, 1-94 (1983).
  * (11) A. Bramon, A. Grau and G. Pancheri, Phys. Lett. B 283, 416-420 (1992).
  * (12) N. Ikeno, M. Bayar and E. Oset, Eur. Phys. J. C 81, 377 (2021).
  * (13) R. L. Workman et al. [Particle Data Group], PTEP 2022, 083C01 (2022).
  * (14) J. X. Lin, J. T. Li, S. J. Jiang, W. H. Liang and E. Oset, Eur. Phys. J. C 81, 1017 (2021).
  * (15) M. Ablikim et al. [BESIII], Phys. Rev. D 95, 032002 (2017).
  * (16) W. H. Liang, J. J. Xie and E. Oset, Eur. Phys. J. C 76, 700 (2016).
  * (17) V. R. Debastiani, W. H. Liang, J. J. Xie and E. Oset, Phys. Lett. B 766, 59-64 (2017).


,

#
The $U(1)_{L_{\mu}-L_{\tau}}$ breaking phase transition, muon $g-2$, dark
matter, collider and gravitational wave

Jie Wang, Jinghong Ma, Jing Gao, Xiao-Fang Han∗00footnotetext: *)
Corresponding author. Email address<EMAIL_ADDRESS>Department of Physics,
Yantai University, Yantai 264005, P. R. China

###### Abstract

Combining the dark matter and muon $g-2$ anomaly, we study the
$U(1)_{L_{\mu}-L_{\tau}}$ breaking phase transition, gravitational wave
spectra, and the direct detection at the LHC in an extra
$U(1)_{L_{\mu}-L_{\tau}}$ gauge symmetry extension of the standard model. The
new fields includes vector-like leptons ($E_{1},~{}E_{2},~{}N$),
$U(1)_{L_{\mu}-L_{\tau}}$ breaking scalar $S$ and gauge boson $Z^{\prime}$, as
well as the dark matter candidate $X_{I}$ and its heavy partner $X_{R}$. A
joint explanation of the dark matter relic density and muon $g-2$ anomaly
excludes the region where both $min(m_{E_{1}},m_{E_{2}},m_{N},m_{X_{R}})$ and
$min(m_{Z^{\prime}},m_{S})$ are much larger than $m_{X_{I}}$. In the parameter
space accommodating the DM relic density and muon $g-2$ anomaly, the model can
achieve a first order $U(1)_{L_{\mu}-L_{\tau}}$ breaking phase transition,
whose strength is sensitive to the parameters of Higgs potential. The
corresponding gravitational wave spectra can reach the sensitivity of
U-DECIGO. In addition, the direct searches at the LHC impose stringent bound
on the mass spectra of the vector-like leptons and dark matter.

###### pacs:

12.60.Fr, 14.80.Ec, 14.80.Bn

##  I introduction

An extra $U(1)_{L_{\mu}-L_{\tau}}$ gauge symmetry extension of the standard
model (SM) is anomaly-free and naturally breaks Lepton Flavour Universality
(LFU) because the $U(1)_{L_{\mu}-L_{\tau}}$ gauge boson couples only to
$\mu(\tau)$ but not to e. The model was originally formulated by He, Joshi,
Lew, and Volkas lu-lt . Thereafter, this type of $U(1)_{L_{\mu}-L_{\tau}}$
model has been modified from its minimal version, and many variants have been
proposed in the context of different phenomenological purpose, such as muon
$g-2$ anomaly Ma:2001md ; Baek:2001kca ; Heeck:2011wj ; Harigaya:2013twa ;
Altmannshofer:2016brv ; Biswas:2016yan ; Biswas:2016yjr ; Banerjee:2018eaf ,
dark matter (DM) puzzle Baek:2008nz ; Das:2013jca ; Patra:2016shz ;
Biswas:2016yan ; Biswas:2016yjr ; Biswas:2017ait ; Foldenauer:2018zrz , and
$b\to s\mu^{+}\mu^{-}$ anomaly bs1 ; bs2 ; bs3 ; bs4 ; bs5 ; bs6 ; bs7 ; bs8 ;
bs9 ; bs10 ; bs11 ; bs12  etc.

In this paper, we will combine the muon $g-2$ anomaly and DM observables, and
examine the $U(1)_{L_{\mu}-L_{\tau}}$ breaking phase transition, gravitational
wave (GW) signatures, and the exclusion of the LHC direct searches in an extra
$U(1)_{L_{\mu}-L_{\tau}}$ gauge symmetry extension of SM. The model was
proposed by one of our authors in Liu:2023sey , in which the new particles
include vector-like leptons ($E_{1},~{}E_{2},~{}N$), $U(1)_{L_{\mu}-L_{\tau}}$
breaking scalar $S$ and gauge boson $Z^{\prime}$, as well as the dark matter
candidate $X_{I}$ and its heavy partner $X_{R}$. When the phase transition is
first order, GW could be generated and detected in current and future GW
experiments, such as LISA lisa , Taiji taiji , TianQin tianqin , Big Bang
Observer (BBO) bbodecigo , DECi-hertz Interferometer GW Observatory (DECIGO)
bbodecigo  and Ultimate-DECIGO (U-DECIGO) udecigo . In addition, the null
results of the LHC direct searches could exclude some parameter space
achieving a first order phase transition (FOPT).

##  II the model

Under the local $U(1)_{L_{\mu}-L_{\tau}}$, the second (third) generation left-
handed lepton doublet and right-handed singlet, $L_{\mu}$, $\mu_{R}$ ,
($L_{\tau}$, $\tau_{R}$), are charged with charge 1 (-1). To obtain the mass
of $U(1)_{L_{\mu}-L_{\tau}}$ gauge boson $Z^{\prime}$, a complex singlet
scalar ${\cal S}$ is required to break the the $U(1)_{L_{\mu}-L_{\tau}}$
symmetry. Another complex singlet scalar $X$ with $U(1)_{L_{\mu}-L_{\tau}}$
charge is introduced whose lighter component may be as a candidate of DM. Also
we add vector-like lepton doublet fields ($E^{"}_{L,R}$) and singlet fields
($E^{\prime}_{L,R}$) which mediate the $X$ interactions to the muon lepton,
and contribute to the muon $g-2$. The quantum numbers of these field under the
gauge group $SU(3)_{C}\times SU(2)_{L}\times U(1)_{Y}\times
U(1)_{L_{\mu}-L_{\tau}}$ are shown in Table 1.

Table 1: The $U(1)_{L_{\mu}-L_{\tau}}$ quantum numbers of the new fields. | SU(3)c | SU(2)L | U(1)Y | U(1)${}_{L_{\mu}-L_{\tau}}$
---|---|---|---|---
$X$ | 1 | 1 | $0$ | $q_{x}$
${\cal S}$ | 1 | 1 | $0$ | $-2q_{x}$
$E^{"}_{L,R}$ | 1 | 2 | $-1/2$ | $1-q_{x}$
$E^{\prime}_{L,R}$ | 1 | 1 | $-1$ | $1-q_{x}$

The Lagrangian is written as

$\displaystyle{\cal L}$ $\displaystyle={\cal L}_{SM}-{1\over
4}Z^{\prime}_{\mu\nu}Z^{\prime\mu\nu}+g_{Z^{\prime}}Z^{\prime\mu}(\bar{\mu}\gamma_{\mu}\mu+\bar{\nu}_{\mu_{L}}\gamma_{\mu}\nu_{\mu_{L}}-\bar{\tau}\gamma_{\mu}\tau-\bar{\nu}_{\tau_{L}}\gamma_{\mu}\nu_{\tau_{L}})$
$\displaystyle+\bar{E"}(i\centernot D)E"+\bar{E^{\prime}}(i\centernot
D)E^{\prime}+(D_{\mu}X^{\dagger})(D^{\mu}X)+(D_{\mu}{\cal
S}^{\dagger})(D^{\mu}{\cal S})$ $\displaystyle-V+\mathcal{L}_{\rm Y}.$ (1)

Where
$Z^{\prime}_{\mu\nu}=\partial_{\mu}Z^{\prime}_{\nu}-\partial_{\nu}Z^{\prime}_{\mu}$
is the field strength tensor, $D_{\mu}$ is the covariant derivative, and
$g_{Z^{\prime}}$ is the gauge coupling constant of the
$U(1)_{L_{\mu}-L_{\tau}}$. $V$ and $\mathcal{L}_{\rm Y}$ denote the scalar
potential and Yukawa interactions.

The scalar potential $V$ containing the SM Higgs parts can be given by

$\displaystyle V$ $\displaystyle=$
$\displaystyle-\mu_{h}^{2}(H^{\dagger}H)-\mu_{S}^{2}({\cal S}^{\dagger}{\cal
S})+m_{X}^{2}(X^{\dagger}X)+\left[\mu X^{2}{\cal S}+\rm h.c.\right]$ (2)
$\displaystyle+\lambda_{H}(H^{\dagger}H)^{2}+\lambda_{S}({\cal
S}^{\dagger}{\cal S})^{2}+\lambda_{X}(X^{\dagger}X)^{2}+\lambda_{SX}({\cal
S}^{\dagger}{\cal S})(X^{\dagger}X)$
$\displaystyle+\lambda_{HS}(H^{\dagger}H)({\cal S}^{\dagger}{\cal
S})+\lambda_{HX}(H^{\dagger}H)(X^{\dagger}X)$

with

$H=\left(\begin{array}[]{c}G^{+}\\\
\frac{1}{\sqrt{2}}\,(h_{1}+v_{h}+iG)\end{array}\right)\,,{\cal
S}={1\over\sqrt{2}}\left(h_{2}+v_{S}+i\omega\right)\,,X={1\over\sqrt{2}}\left(X_{R}+iX_{I}\right)\,.$
(3)

Where $v_{h}=246$ GeV and $v_{S}$ are respectively vacuum expectation values
(VeVs) of $H$ and ${\cal S}$, and the $X$ field has no VeV. One can determine
the mass parameters $\mu^{2}_{h}$ and $\mu^{2}_{S}$ of Eq. (2) by the
potential minimization conditions,

$\begin{split}&\quad\mu_{h}^{2}=\lambda_{H}v_{h}^{2}+{1\over
2}\lambda_{HS}v_{S}^{2},\\\ &\quad\mu_{S}^{2}=\lambda_{S}v_{S}^{2}+{1\over
2}\lambda_{HS}v_{h}^{2}.\\\ \end{split}$ (4)

After ${\cal S}$ acquires the VeV, the $\mu$ term makes the complex scalar $X$
split into two real scalar fields ($X_{R}$, $X_{I}$), and their masses are

$\displaystyle m_{X_{R}}^{2}=m_{X}^{2}+{1\over 2}\lambda_{HX}v_{H}^{2}+{1\over
2}\lambda_{SX}v_{S}^{2}+\sqrt{2}\mu v_{S}$ $\displaystyle
m_{X_{I}}^{2}=m_{X}^{2}+{1\over 2}\lambda_{HX}v_{H}^{2}+{1\over
2}\lambda_{SX}v_{S}^{2}-\sqrt{2}\mu v_{S}.$ (5)

After the $U(1)_{L_{\mu}-L_{\tau}}$ is broken, there is still remnant $Z_{2}$
symmetry, which guarantee the lightest component $X_{I}$ to be as a candidate
of DM.

The $\lambda_{HS}$ and $\lambda_{HX}$ terms will lead to the couplings of 125
GeV Higgs ($h$) and DM. To suppress the stringent constraints from the DM
direct detection and indirect detection experiments, we simply assume that the
$hX_{I}X_{I}$ coupling is absent, namely taking $\lambda_{HX}=0$ and
$\lambda_{HS}=0$. Thus, the 125 GeV Higgs $h$ is purely from $h_{1}$ and extra
CP-even Higgs $S$ is purely from $h_{2}$. Their masses are given by

$m_{h}^{2}=2\lambda_{H}v_{h}^{2},~{}~{}~{}m_{S}^{2}=2\lambda_{S}v_{S}^{2}.$
(6)

After ${\cal S}$ acquires VeV, the $U(1)_{L_{\mu}-L_{\tau}}$ gauge boson
$Z^{\prime}$ obtains a mass,

$m_{Z}^{\prime}=2g_{Z}^{\prime}\mid q_{x}\mid v_{S}.$ (7)

The Yukawa interactions $\mathcal{L}_{\rm Y}$ can be written as

$\displaystyle-\mathcal{L}_{\rm Y,mass}=$ $\displaystyle
m_{1}\overline{E^{\prime}_{L}}E^{\prime}_{R}+m_{2}\overline{E^{"}_{L}}E^{"}_{R}+\kappa_{1}\overline{\mu_{R}}XE^{\prime}_{L}+\kappa_{2}\overline{L_{\mu}}XE^{"}_{R}$
$\displaystyle+\sqrt{2}y_{1}\overline{E^{"}_{L}}HE^{\prime}_{R}+\sqrt{2}y_{2}\overline{E^{"}_{R}}HE^{\prime}_{L}+\frac{\sqrt{2}m_{\mu}}{v}\overline{L_{\mu}}H\mu_{R}+{\rm
h.c.}.$ (8)

After the Electroweak symmetry breaking, the vector-like lepton masses are
given by

$\displaystyle M_{E}=\begin{pmatrix}m_{1}&y_{2}v_{h}\\\
y_{1}v_{h}&m_{2}\end{pmatrix}.$ (9)

By making a bi-unitary transformation with the rotation matrices for the
right-handed fields and the left-handed fields,

$\displaystyle U_{L}=\begin{pmatrix}c_{L}&-s_{L}\\\
s_{L}&c_{L}\end{pmatrix},\quad U_{R}=\begin{pmatrix}c_{R}&-s_{R}\\\
s_{R}&c_{R}\end{pmatrix},$ (10)

where $c_{L,R}^{2}+s_{L,R}^{2}=1$, we can diagonalize the mass matrix for the
vector-like lepton,

$\displaystyle
U_{L}^{\dagger}M_{E}U_{R}=\mathrm{diag}\left(m_{E_{1}},m_{E_{2}}\right).$ (11)

The $E_{1}$ and $E_{2}$ are the mass eigenstates of charged vector-like
leptons, and the mass of neutral vector-like lepton $N$ is

$m_{N}=m_{2}=m_{E_{2}}c_{L}c_{R}+m_{E_{1}}s_{L}s_{R}.$ (12)

The $E_{1}$ and $E_{2}$ can mediate $X_{R}$ and $X_{I}$ interactions to muon
lepton,

$-\mathcal{L}_{\rm
X}\supset\frac{1}{\sqrt{2}}(X_{R}+iX_{I})\left[\bar{\mu}_{R}(\kappa_{1}c_{L}E_{1L}-\kappa_{1}s_{L}E_{2L})+\bar{\mu}_{L}(\kappa_{2}s_{R}E_{1R}+\kappa_{2}c_{R}E_{2R})\right]+h.c.~{},$
(13)

and have the couplings to the 125 GeV Higgs,

$\displaystyle-\mathcal{L}_{\rm h}$ $\displaystyle\supset$
$\displaystyle\frac{m_{E_{1}}(c_{L}^{2}s_{R}^{2}+c_{R}^{2}s_{L}^{2})-2m_{E_{2}}s_{L}c_{L}s_{R}c_{R}}{v_{h}}~{}h\bar{E}_{1}E_{1},$
(14)
$\displaystyle+\frac{m_{E_{2}}(s_{L}^{2}c_{R}^{2}+c_{L}^{2}s_{R}^{2})-2m_{E_{1}}s_{L}c_{L}s_{R}c_{R}}{v_{h}}~{}h\bar{E}_{2}E_{2}.$

##  III Dark matter and muon $g-2$

We fix on $m_{h}=$ 125 GeV, $v_{h}=$ 246 GeV, $q_{x}=-1$, $\lambda_{HS}=0$,
and $\lambda_{HX}=0$, and take $g_{Z^{\prime}}$, $m_{Z^{\prime}}$,
$\lambda_{X}$, $\lambda_{SX}$, $m_{S}$, $m_{X_{R}}$, $m_{X_{I}}$, $m_{E_{1}}$,
$m_{E_{2}}$, $s_{L}$, $s_{R}$, $\kappa_{1}$, and $\kappa_{2}$ as the free
parameters. To maintain the perturbativity, we conservatively choose

$\displaystyle\mid\lambda_{SX}\mid\leq 4\pi,~{}~{}\mid\lambda_{X}\mid\leq
4\pi,$
$\displaystyle-\frac{1}{2}\leq\kappa_{1}\leq\frac{1}{2},~{}~{}-\frac{1}{2}\leq\kappa_{2}\leq\frac{1}{2},$
(15)

and take the mixing parameters $s_{L}$ and $s_{R}$ as

$-\frac{1}{\sqrt{2}}\leq
s_{L}\leq\frac{1}{\sqrt{2}},~{}~{}-\frac{1}{\sqrt{2}}\leq
s_{R}\leq\frac{1}{\sqrt{2}}.$ (16)

The mass parameters are scanned over in the following ranges:

$\displaystyle 60~{}{\rm GeV}\leq m_{X_{I}}\leq 500~{}{\rm
GeV},~{}~{}~{}m_{X_{I}}\leq m_{X_{R}}\leq 500~{}{\rm GeV},$ $\displaystyle
m_{X_{I}}\leq m_{E_{1}}\leq 500~{}{\rm GeV},~{}~{}~{}m_{X_{I}}\leq
m_{E_{2}}\leq 500~{}{\rm GeV},$ $\displaystyle 100~{}{\rm GeV}\leq
m_{Z^{\prime}}\leq 500~{}{\rm GeV},~{}~{}~{}100~{}{\rm GeV}\leq m_{S}\leq
500~{}{\rm GeV}.$ (17)

The mass of neutral vector-like lepton $N$ is a function of $m_{E_{1}}$,
$m_{E_{2}}$, $s_{L}$ and $s_{R}$, and $m_{X_{I}}<m_{N}$ is imposed. We require
0 $<g_{Z^{\prime}}/m_{Z^{\prime}}\leq$ (550 GeV)-1 to be consistent with the
bound of the neutrino trident process trident .

The potential stability in Eq. (2) requires the following condition,

$\displaystyle\lambda_{H}\geq 0\,,\quad\lambda_{S}\geq
0\,,\quad\lambda_{X}\geq 0\,,\quad$
$\displaystyle\lambda_{HS}\geq-2\sqrt{\lambda_{H}\,\lambda_{S}}\,,\quad\lambda_{HX}\geq-2\sqrt{\lambda_{H}\,\lambda_{X}}\,,\quad\lambda_{SX}\geq-2\sqrt{\lambda_{S}\,\lambda_{X}}\,,\quad$
$\displaystyle\sqrt{\lambda_{HS}+2\sqrt{\lambda_{H}\,\lambda_{S}}}~{}\sqrt{\lambda_{HX}+2\sqrt{\lambda_{H}\,\lambda_{X}}}~{}\sqrt{\lambda_{SX}+2\sqrt{\lambda_{S}\,\lambda_{X}}}$
$\displaystyle+2\,\sqrt{\lambda_{H}\lambda_{S}\lambda_{X}}+\lambda_{HS}\sqrt{\lambda_{X}}+\lambda_{HX}\sqrt{\lambda_{S}}+\lambda_{SX}\sqrt{\lambda_{H}}\geq
0\,.$ (18)

The one-loop diagrams with the vector-like lepton can give additional
corrections to the oblique parameters ($S,~{}T,~{}U$), which can be calculated
as in Refs. stu ; mc-stu ; 1305.4712 . Taking the recent fit results of Ref.
pdg2020 , we use the following values of $S$, $T$, and $U$,

$S=-0.01\pm 0.10,~{}~{}T=0.03\pm 0.12,~{}~{}U=0.02\pm 0.11,$ (19)

with the correlation coefficients

$\rho_{ST}=0.92,~{}~{}\rho_{SU}=-0.80,~{}~{}\rho_{TU}=-0.93.$ (20)

Also the one-loop diagrams of the charged vector-like leptons $E_{1}$ and
$E_{2}$ can contribute to the $h\to\gamma\gamma$ decay, and the bound of the
diphoton signal strength of the 125 GeV Higgs is imposed pdg2020 ,

$\mu_{\gamma\gamma}=1.11^{+0.10}_{-0.09}.$ (21)

In the model, the dominant corrections to the muon $g-2$ are from the one-loop
diagrams with the vector-like leptons ($E_{1},~{}E_{2})$ and scalar fields
($X_{R}$ and $X_{I}$), which are approximately calculated as in Refs.
1305.3522 ; 1906.11297 ; 0902.3360

$\displaystyle\Delta a_{\mu}$ $\displaystyle=$
$\displaystyle\frac{1}{32\pi^{2}}m_{\mu}\left(\kappa_{1}c_{L}\kappa_{2}s_{R}H(m_{E_{1}},m_{X_{R}})-\kappa_{1}s_{L}\kappa_{2}c_{R}H(m_{E_{2}},m_{X_{R}})\right.$
(22)
$\displaystyle\left.+\kappa_{1}c_{L}\kappa_{2}s_{R}H(m_{E_{1}},m_{X_{I}})-\kappa_{1}s_{L}\kappa_{2}c_{R}H(m_{E_{2}},m_{X_{I}})\right),$

where the function

$H(m_{f},m_{\phi})=\frac{m_{f}}{m_{\phi}^{2}}\frac{(r^{2}-4r+2\log{r}+3)}{(r-1)^{3}}$
(23)

with $r=\frac{m_{f}^{2}}{m_{\phi}^{2}}$. The combined average for the muon
$g-2$ with Fermilab E989 fermig2  and Brookhaven E821 E821 , the difference
from the SM prediction becomes

$\displaystyle\Delta a_{\mu}=a_{\mu}^{exp}-a_{\mu}^{SM}=(25.1\pm 5.9)\times
10^{-10},$ (24)

which shows $4.2\sigma$ discrepancy from the SM. Very recently, on August 10,
2023, the E989 experiment at Fermilab released an update regarding the
measurement from Run-2 and Run-3 mug2-2023 . The new combined value yields a
deviation of

$\displaystyle\Delta a_{\mu}=a_{\mu}^{exp}-a_{\mu}^{SM}=(24.9\pm 4.8)\times
10^{-10},$ (25)

which leads to a $5.1\sigma$ discrepancy. Whereas the recent lattice
calculation mug2-lat  and the experiment determination mug2-had  of the hadron
vacuum polarization contribution to the muon $g-2$ point the value closer to
the SM prediction, and hence the tension relaxes to a few sigma level.

If kinematically allowed, the DM pair-annihilation processes includes
$X_{I}X_{I}\to\mu^{+}\mu^{-},~{}Z^{\prime}Z^{\prime},~{}SS$. In addition, a
small mass splitting between the DM and the other new particles
($E_{1},~{}E_{2},~{}N$, $X_{R}$) can lead to coannihilation. The Planck
collaboration reported the relic density of cold DM in the universe,
$\Omega_{c}h^{2}=0.1198\pm 0.0015$ planck , and the theoretical prediction of
the model is calculated by micrOMEGAs-5.2.13 micomega .

Figure 1: The surviving samples explaining the DM relic density and muon $g-2$
anomaly while satisfying the constraints from theory, oblique parameter, and
125 GeV Higgs diphoton signal.

After imposing the constraints of theory, oblique parameters, and the diphoton
signal data of the 125 GeV Higgs, we project the samples accommodating the DM
relic density and muon $g-2$ anomaly within $2\sigma$ ranges in Fig. 1. From
Fig. 1 we find that the correct DM relic density can be obtained for most of
the parameter region of
$min(m_{E_{1}},~{}m_{E_{2}},~{}m_{N},~{}m_{X_{R}})-m_{X_{I}}<350$ GeV and -400
GeV $<min(m_{Z^{\prime}},m_{S})-m_{X_{I}}<$ 400 GeV. The DM relic density is
mainly produced via the DM pair-annihilation processes
$X_{I}X_{I}\to~{}Z^{\prime}Z^{\prime},~{}SS$ for the region of
$min(m_{Z^{\prime}},m_{S})<m_{X_{I}}$, the coannihilation processes for the
region of $min(m_{E_{1}},~{}m_{E_{2}},~{}m_{N},~{}m_{X_{R}})$ close to
$m_{X_{I}}$, and $X_{I}X_{I}\to\mu^{+}\mu^{-}$ for the region of both
$min(m_{E_{1}},~{}m_{E_{2}},~{}m_{N},~{}m_{X_{R}})>m_{X_{I}}$ and
$min(m_{Z^{\prime}},m_{S})>m_{X_{I}}$. However, once the explanation of muon
$g-2$ anomaly is simultaneously required, most of region of
$min(m_{E_{1}},~{}m_{E_{2}},~{}m_{N},~{}m_{X_{R}})>m_{X_{I}}$ and
$min(m_{Z^{\prime}},m_{S})>m_{X_{I}}$ is ruled out. This is because the muon
$g-2$ anomaly favors small interactions between the vector-like leptons and
muon mediated by $X_{I}$, which leads to $X_{I}X_{I}\to\mu^{+}\mu^{-}$ process
to fail to produce the correct DM relic density.

The $X_{I}$ does not couple to the SM quark, and its couplings to the muon
lepton and vector-like leptons are restricted by the muon $g-2$ anomaly.
Therefore, the model can accommodate the bound from the DM direct detection
naturally.

##  IV $U(1)_{L_{\mu}-L_{\tau}}$ breaking phase transition

At high temperatures, the global minimum of the finite-temperature effective
potential is at the origin, i.e. $SU(2)_{L}\times U(1)_{Y}\times
U(1)_{L_{\mu}-L_{\tau}}$ is unbroken. When the temperature drops, the
potential changes and at some point develops a minimum at non-vanishing field
values. The phase transition between the unbroken and the broken phase can
proceed in basically two different ways. In a FOPT, at the critical
temperature $T_{C}$, the two degenerate minima will be at different points in
field space, typically with a potential barrier in between. For a second order
(cross-over) transition, the broken and symmetric minimum are not degenerate
until they are at the same point in field space. In this paper we focus on a
first-order $U(1)_{L_{\mu}-L_{\tau}}$ breaking phase transition.

###  IV.1 The thermal effective potential

In order to examine phase transition, we first take $h_{1}$, $h_{2}$, and
$X_{r}$ as the field configurations, and obtain the field dependent masses of
the scalars ($h$, $S$, $X_{R}$, $X_{I}$), the Goldstone boson
($G,~{}\omega,~{}G^{\pm}$), the gauge boson, and fermions. The field dependent
masses of scalars are given

$\displaystyle{\hat{m}}^{2}_{h,S,X_{R}}$
$\displaystyle=\rm{eigenvalues}(\widehat{\mathcal{M}^{2}_{P}})\ ,$ (26)
$\displaystyle{\hat{m}}^{2}_{G,\omega,X_{I}}$
$\displaystyle=\rm{eigenvalues}(\widehat{\mathcal{M}^{2}_{A}})\ ,$ (27)
$\displaystyle{\hat{m}}^{2}_{G^{\pm}}$
$\displaystyle=\lambda_{H}h_{1}^{2}-\lambda_{H}v_{h}^{2}\ ,$ (28)

with

$\displaystyle\widehat{\mathcal{M}^{2}_{P}}_{11}$
$\displaystyle=3\lambda_{H}h^{2}_{1}-\lambda_{H}v_{h}^{2}$
$\displaystyle\widehat{\mathcal{M}^{2}_{P}}_{22}$
$\displaystyle=-\lambda_{S}v_{S}^{2}+3\lambda_{S}h_{2}^{2}+{\lambda_{SX}\over
2}X_{r}^{2}$ $\displaystyle\widehat{\mathcal{M}^{2}_{P}}_{33}$
$\displaystyle=m_{X}^{2}+\sqrt{2}\mu h_{2}+{\lambda_{SX}\over
2}h_{2}^{2}+3\lambda_{X}X_{r}^{2}$
$\displaystyle\widehat{\mathcal{M}^{2}_{P}}_{23}$
$\displaystyle=\widehat{\mathcal{M}^{2}_{P}}_{32}=\sqrt{2}\mu
X_{r}+\lambda_{SX}h_{2}X_{r}$
$\displaystyle\widehat{\mathcal{M}^{2}_{P}}_{12}$
$\displaystyle=\widehat{\mathcal{M}^{2}_{P}}_{21}=\widehat{\mathcal{M}^{2}_{P}}_{13}=\widehat{\mathcal{M}^{2}_{P}}_{31}=0$
$\displaystyle\widehat{\mathcal{M}^{2}_{A}}_{11}$
$\displaystyle=\lambda_{H}h^{2}_{1}-\lambda_{H}v_{h}^{2}$
$\displaystyle\widehat{\mathcal{M}^{2}_{A}}_{22}$
$\displaystyle=-\lambda_{S}v_{S}^{2}+\lambda_{S}h^{2}_{2}+{\lambda_{SX}\over
2}X_{r}^{2}$ $\displaystyle\widehat{\mathcal{M}^{2}_{A}}_{33}$
$\displaystyle=m_{X}^{2}-\sqrt{2}\mu
h_{2}+\lambda_{X}X^{2}_{r}+{\lambda_{SX}\over 2}h_{2}^{2}$
$\displaystyle\widehat{\mathcal{M}^{2}_{A}}_{23}$
$\displaystyle=\widehat{\mathcal{M}^{2}_{A}}_{32}=-\sqrt{2}\mu X_{r}$
$\displaystyle\widehat{\mathcal{M}^{2}_{A}}_{12}$
$\displaystyle=\widehat{\mathcal{M}^{2}_{A}}_{21}=\widehat{\mathcal{M}^{2}_{A}}_{13}=\widehat{\mathcal{M}^{2}_{A}}_{31}=0.$
(29)

The field dependent masses of gauge boson are

$\displaystyle{\hat{m}}^{2}_{W^{\pm}}={1\over 4}g^{2}h^{2}_{1},$
$\displaystyle{\hat{m}}^{2}_{Z}={1\over 4}(g^{2}+g^{\prime 2})h^{2}_{1},$ (30)
$\displaystyle\quad{\hat{m}}^{2}_{\gamma}=0,$
$\displaystyle{\hat{m}}^{2}_{Z^{\prime}}=q_{x}^{2}g_{Z^{\prime}}^{2}(4h^{2}_{2}+X^{2}_{r}),$
(31)

The field dependent masses of vector-like lepton are

$\displaystyle\widehat{\mathcal{M}}^{2}_{E_{1},E_{2},\mu}$
$\displaystyle=\rm{eigenvalues}(\widehat{\mathcal{M}_{E}}\widehat{\mathcal{M}_{E}}^{T})$
(32)

with

$\displaystyle\widehat{\mathcal{M}_{E}}=\begin{pmatrix}m_{1}&\quad
y_{2}h_{1}&\quad\kappa_{1}X_{r}\\\\[5.0pt] y_{1}h_{1}&\quad m_{2}&\quad
0\\\\[5.0pt]
0&\quad\kappa_{2}X_{r}&\quad\frac{m_{\mu}h_{1}}{v_{h}}\end{pmatrix}.$ (33)

For the quarks of SM, we only consider the top quark,

$\displaystyle\widehat{\mathcal{M}^{2}_{t}}$
$\displaystyle=y_{t}^{2}h_{1}^{2}$ (34)

with $y_{t}=m_{t}/v_{h}$.

In order to examine the $U(1)_{L_{\mu}-L_{\tau}}$ breaking phase transition,
we need to study the thermal effective potential $V_{eff}$ in terms of the
classical fields ($h_{1},h_{2},X_{r}$), which is composed of four parts:

$\displaystyle V_{eff}(h_{1},h_{2},X_{r},T)=$ $\displaystyle
V_{0}(h_{1},h_{2},X_{r})+V_{CW}(h_{1},h_{2},X_{r})+V_{CT}(h_{1},h_{2},X_{r})$
$\displaystyle+V_{T}(h_{1},h_{2},X_{r},T)+V_{ring}(h_{1},h_{2},X_{r},T).$ (35)

$V_{0}$ is the tree-level potential, $V_{CW}$ is the Coleman-Weinberg (CW)
potential vcw , $V_{CT}$ is the counter term, $V_{T}$ is the thermal
correction vloop , and $V_{ring}$ is the resummed daisy corrections vring1 ;
vring2 . In this paper, we calculate $V_{eff}$ in the Landau gauge.

The tree-level potential $V_{0}$ in terms of their classical fields
($h_{1},h_{2},X_{r}$) from the Eq. (2),

$\displaystyle V_{0}(h_{1},h_{2},X_{r})$ $\displaystyle=$
$\displaystyle-\frac{\mu_{h}^{2}}{2}h_{1}^{2}-\frac{\mu_{S}^{2}}{2}h_{2}^{2}+\frac{m_{X}^{2}}{2}X_{r}^{2}+\frac{\mu}{\sqrt{2}}h_{2}X_{r}^{2}$
(36)
$\displaystyle+\frac{\lambda_{H}}{4}h_{1}^{4}+\frac{\lambda_{S}}{4}h_{2}^{4}+\frac{\lambda_{X}}{4}X_{r}^{4}+\frac{\lambda_{SX}}{4}X_{r}^{2}h_{2}^{2}.$

The CW potential in the $\overline{\rm MS}$ scheme at 1-loop level has the
form vcw :

$V_{\rm
CW}(h_{1},h_{2},X_{r})=\sum_{i}(-1)^{2s_{i}}n_{i}\frac{{\hat{m}}_{i}^{4}(h_{1},h_{2},X_{r})}{64\pi^{2}}\left[\ln\frac{{\hat{m}}_{i}^{2}(h_{1},h_{2},X_{r})}{Q^{2}}-C_{i}\right]\;\
,$ (37)

where
$i=h,S,X_{R},X_{I},G,\omega,G^{\pm},W^{\pm},Z,Z^{\prime},t,E_{1},E_{2},\mu$,
and $s_{i}$ is the spin of particle i. $Q$ is a renormalization scale, and we
take $Q=m_{S}$. The constants $C_{i}=\frac{3}{2}$ for scalars or fermions and
$C_{i}=\frac{5}{6}$ for gauge bosons. $n_{i}$ is the number of degree of
freedom,

$\displaystyle n_{h}=n_{S}=n_{X_{R}}=n_{X_{I}}=n_{G}=n_{\omega}=1$
$\displaystyle n_{G^{\pm}}=2,~{}~{}n_{W^{\pm}}=6,~{}n_{Z}=n_{Z^{\prime}}=3$
$\displaystyle n_{t}=12,~{}~{}n_{E_{1}}=n_{E_{2}}=n_{\mu}=4.$ (38)

With $V_{CW}$ being included in the potential, the minimization conditions of
scalar potential and the CP-even mass matrix will be shifted slightly. To
maintain these relations, the counter terms $V_{ct}$ should be added,

$\displaystyle V_{\rm CT}$ $\displaystyle=\delta m_{1}^{2}h_{1}^{2}+\delta
m_{2}^{2}h_{2}^{2}+\delta
m_{X}^{2}X_{r}^{2}+\delta\lambda_{H}h_{1}^{4}+\delta\lambda_{S}h_{2}^{4}+\delta\lambda_{X}X_{r}^{4}$
$\displaystyle+\delta\mu X_{r}^{2}h_{2}+\delta\lambda_{SX}h_{2}^{2}X_{r}^{2}.$
(39)

The relevant coefficients are determined by

$\frac{\partial V_{\rm CT}}{\partial h_{1}}=-\frac{\partial V_{\rm
CW}}{\partial h_{1}}\;,\quad\frac{\partial V_{\rm CT}}{\partial
h_{2}}=-\frac{\partial V_{\rm CW}}{\partial h_{2}},\;\quad\frac{\partial
V_{\rm CT}}{\partial X_{r}}=-\frac{\partial V_{\rm CW}}{\partial X_{r}},$ (40)
$\displaystyle\frac{\partial^{2}V_{\rm CT}}{\partial h_{1}\partial h_{1}}$
$\displaystyle=-\frac{\partial^{2}V_{\rm CW}}{\partial h_{1}\partial
h_{1}}\;,\quad\frac{\partial^{2}V_{\rm CT}}{\partial h_{2}\partial
h_{2}}=-\frac{\partial^{2}V_{\rm CW}}{\partial h_{2}\partial
h_{2}}\;,\quad\frac{\partial^{2}V_{\rm CT}}{\partial X_{r}\partial
X_{r}}=-\frac{\partial^{2}V_{\rm CW}}{\partial X_{r}\partial X_{r}}\;,$
$\displaystyle\frac{\partial^{2}V_{\rm CT}}{\partial h_{1}\partial h_{2}}$
$\displaystyle=-\frac{\partial^{2}V_{\rm CW}}{\partial h_{1}\partial
h_{2}}\;,\quad\frac{\partial^{2}V_{\rm CT}}{\partial h_{1}\partial
X_{r}}=-\frac{\partial^{2}V_{\rm CW}}{\partial h_{1}\partial
X_{r}}\;,\quad\frac{\partial^{2}V_{\rm CT}}{\partial h_{2}\partial
X_{r}}=-\frac{\partial^{2}V_{\rm CW}}{\partial h_{2}\partial X_{r}}\;,\quad$
(41)

which are evaluated at the EW minimum of
$\\{h_{1}=v_{h},h_{2}=v_{S},X_{r}=0\\}$ on both sides. As a result, the VeVs
of $h_{1}$, $h_{2}$, $X_{r}$ and the CP-even mass matrix will not be shifted.
We check that the following relations hold true

$\displaystyle\frac{\partial^{2}V_{\rm CW}}{\partial h_{1}\partial
h_{2}}=0\;,\quad\frac{\partial^{2}V_{\rm CW}}{\partial h_{1}\partial
X_{r}}=0\;.$ (42)

For the left seven equations, there are eight parameters to be fixed, so that
one renormalization constant is left for determination. We choose to use
$\delta m_{1}^{2}$, $\delta m_{2}^{2}$, $\delta m_{X}^{2}$,
$\delta\lambda_{H}$, $\delta\lambda_{S}$, $\delta\lambda_{X}$,
$\delta\lambda_{SX}$, and set $\delta\mu=0$. It is a well-known problem that
the second derivative of the CW potential in the vacuum suffers from
logarithmic divergences originating from the vanishing Goldstone masses. In
order to avoid the problems with infrared divergent Goldstone contributions,
we simply remove the Goldstone corrections in the renormalization conditions
following the approaches of 1107.3559 ; Espinosa:2011ax . This is simply a
change of renormalization conditions and the shift it causes in the potential
shape is negligible.

The thermal contributions $V_{T}$ to the potential can be written as vloop

$V_{\rm
th}(h_{1},h_{2},X_{r},T)=\frac{T^{4}}{2\pi^{2}}\,\sum_{i}n_{i}J_{B,F}\left(\frac{{\hat{m}}_{i}^{2}(h_{1},h_{2},X_{r})}{T^{2}}\right)\;,$
(43)

where
$i=h,S,X_{R},X_{I},G,\omega,G^{\pm},W^{\pm},Z,Z^{\prime},t,E_{1},E_{2},\mu$,
and the functions $J_{B,F}$ are

$J_{B,F}(y)=\pm\int_{0}^{\infty}\,dx\,x^{2}\,\ln\left[1\mp{\rm
exp}\left(-\sqrt{x^{2}+y}\right)\right].$ (44)

Finally, the thermal corrections with resumed ring diagrams are given vring1 ;
vring2 .

$V_{\rm
ring}\left(h_{1},h_{2},X_{r},T\right)=-\frac{T}{12\pi}\sum_{i}n_{i}\left[\left(\bar{M}_{i}^{2}\left(h_{1},h_{2},X_{r},T\right)\right)^{\frac{3}{2}}-\left({\hat{m}}_{i}^{2}\left(h_{1},h_{2},X_{r},T\right)\right)^{\frac{3}{2}}\right],$
(45)

where
$i=h,S,X_{R},X_{I},G,\omega,G^{\pm},W^{\pm}_{L},Z_{L},Z^{\prime}_{L},\gamma_{L}$.
The $W^{\pm}_{L},~{}Z_{L},~{}Z^{\prime}_{L}$, and $\gamma_{L}$ are the
longitudinal gauge bosons with
$n_{W^{\pm}_{L}}=2,~{}n_{Z_{L}}=n_{Z^{\prime}_{L}}=n_{\gamma_{L}}=1$. The
thermal Debye masses $\bar{M}_{i}^{2}\left(h_{1},h_{2},X_{r},T\right)$ for the
CP-even and CP-odd scalar fields are the eigenvalues of the full mass matrix,

$\bar{M}_{i}^{2}\left(h_{1},h_{2},X_{r},T\right)={\rm
eigenvalues}\left[\widehat{\mathcal{M}_{X}^{2}}\left(h_{1},h_{2},X_{r}\right)+\Pi_{X}(T)\right],$
(46)

where $\Pi_{X}$ ($X=P,A$) are given by

$\displaystyle(\Pi_{P,A})_{11}$ $\displaystyle=\left[{9g^{2}\over
2}+{3g^{\prime 2}\over
2}+12y_{t}^{2}+4(y_{2}^{2}+y_{1}^{2})+12\lambda_{H}\right]{T^{2}\over 24}$
$\displaystyle(\Pi_{P,A})_{22}$
$\displaystyle=\left[8\lambda_{S}+2\lambda_{SX}+24q_{x}^{2}g^{2}_{Z^{\prime}}\right]{T^{2}\over
24}$ $\displaystyle(\Pi_{P,A})_{33}$
$\displaystyle=\left[8\lambda_{X}+2\lambda_{SX}+6q_{x}^{2}g^{2}_{Z^{\prime}}+2\kappa_{1}^{2}+4\kappa_{2}^{2}\right]{T^{2}\over
24}$ $\displaystyle(\Pi_{P,A})_{12}$
$\displaystyle=(\Pi_{P,A})_{21}=(\Pi_{P,A})_{13}=(\Pi_{P,A})_{31}=(\Pi_{P,A})_{32}=(\Pi_{P,A})_{23}=0.$
(47)

The thermal Debye mass of $G^{\pm}$ is

$\displaystyle\bar{M}_{G^{+}G^{-}}^{2}\left(h_{1},h_{2},X_{r},T\right)=\lambda_{H}h_{1}^{2}-\lambda_{H}v_{h}^{2}+\Pi_{P11}.$
(48)

The Debye masses of longitudinal gauge bosons are

$\displaystyle\bar{M}_{W^{\pm}_{L}}^{2}\left(h_{1},h_{2},X_{r},T\right)$
$\displaystyle={1\over 4}g^{2}h^{2}_{1}+{11\over 6}g^{2}T^{2},$
$\displaystyle\bar{M}_{Z_{L},\gamma_{L}}^{2}\left(h_{1},h_{2},X_{r},T\right)$
$\displaystyle=\frac{1}{8}\left(g^{2}+g^{\prime
2}\right)\left(h_{1}^{2}+\frac{22}{3}T^{2}\right)\pm\frac{1}{2}\Delta,$
$\displaystyle\bar{M}_{Z^{\prime}_{L}}^{2}\left(h_{1},h_{2},X_{r},T\right)$
$\displaystyle=q_{x}^{2}g_{Z^{\prime}}^{2}(4h^{2}_{2}+X^{2}_{r})+{1\over
6}g^{2}_{Z^{\prime}}T^{2}\left[12-12q_{x}+16q_{x}^{2}\right],$ (49)

with $\Delta=\sqrt{\left[\frac{1}{4}\left(g^{2}-g^{\prime
2}\right)\left(h_{1}^{2}+\frac{22}{3}T^{2}\right)\right]^{2}+\frac{1}{4}g^{2}g^{\prime
2}h_{1}^{2}}$.

Figure 2: The surviving samples achieving FOFP and accommodating the DM relic
density, muon $g-2$ anomaly, and various constraints mentioned previously.
$\xi_{C}=\frac{<S>}{T_{C}}$ denotes the strength of the
$U(1)_{L_{\mu}-L_{\tau}}$ breaking FOPT at the critical temperature $T_{C}$.

In Fig. 2, we display the scatter plots achieving FOFP and accommodating the
DM relic density, muon $g-2$ anomaly, and various constraints mentioned
previously. We observe that the strength of FOPT is sensitive to the parameter
$\lambda_{S}$. As $\lambda_{S}$ decreases, the the critical temperature
$T_{C}$ tends to decrease, and the strength of FOPT tends to increase.

##  V Collider and gravitational wave signature

###  V.1 Limits from the collider experiments

Since the $U(1)_{L_{\mu}-L_{\tau}}$ gauge boson $Z^{\prime}$ does not couple
to quarks, the $Z^{\prime}$ is rather hard to produce at the LHC. The vector-
like leptons are mainly pair produced at the LHC via the electroweak processes
mediated by the SM gauge bosons, and the $E_{1,2}\to\mu X_{I}$ and
$N\to\nu_{\mu}X_{I}$ are main decay modes of vector-like leptons. Thus, we
employ the ATLAS analysis of $2\ell+E_{T}^{miss}$ with 139 fb-1 integrated
luminosity data to restrict our model 1908.08215 , which is implemented in the
MadAnalysis5 ma5-1 ; ma5-2 ; ma5-3  with assuming 95% confidence level for the
exclusion limit. The simulations for the samples are performed by
MG5_aMC-3.3.2 Alwall:2014hca  with PYTHIA8 Torrielli:2010aw  and Delphes-3.2.0
deFavereau:2013fsa .

Figure 3: All the samples achieve FOFP and accommodate the DM relic density,
muon $g-2$ anomaly, and various constraints mentioned previously. The circles
and squares are excluded and allowed by the direct searches for
$2\ell+E^{miss}_{T}$ at the LHC.

In Fig. 3, we employ the ATLAS analysis of $2\ell+E_{T}^{miss}$ at the LHC to
restrict the parameter space which has been satisfied by the DM relic density,
muon $g-2$, FOPT, and various constraints discussed above. Fig. 3 indicates
that the DM mass is allowed to be as low as 100 GeV when the value of
$min(m_{E_{1}},~{}m_{E_{2}})-m_{X_{I}}$ is small. This is because the muon
from the vector-like lepton decay has soft energy, and its detection
efficiency is decreased at the LHC. As the DM mass increase, the value of
$min(m_{E_{1}},~{}m_{E_{2}})-m_{X_{I}}$ increases, and the lightest charged
vector-like lepton is allowed to have a more large mass.

###  V.2 Gravitational wave signature

Stochastic GWs are produced during a FOPT via bubble collision, sound waves in
the plasma and the magneto-hydrodynamics turbulence. Since most of the phase
transition energy is pumped into the surrounding fluid shells, making sound
waves the dominant contribution to GWs, we will focus on the GW spectrum from
the sound waves in the plasma.

The sound wave spectra can be expressed as functions of two FOPT parameters
$\beta$ and $\alpha$,

$\displaystyle\frac{\beta}{H_{n}}=T\frac{d(S_{3}(T)/T)}{dT}|_{T=T_{n}}\;,~{}~{}~{}~{}\alpha=\frac{\Delta\rho}{\rho_{R}}=\frac{\Delta\rho}{\pi^{2}g_{\ast}T_{n}^{4}/30}\;.$
(50)

Where $H_{n}$ is the Hubble parameter at the nucleation temperature $T_{n}$,
and $g_{\ast}$ is the effective number of relativistic degrees of freedom.
$\beta$ characterizes roughly the inverse time duration of the strong first-
order PT, and $\alpha$ is defined as the vacuum energy released from the phase
transition normalized by the total radiation energy density $\rho_{R}$ at
$T_{n}$.

The GW spectrum from the sound waves can be expressed by gw-sw

$\displaystyle\Omega_{\textrm{sw}}h^{2}$ $\displaystyle\ =\ $ $\displaystyle
2.65\times
10^{-6}\left(\frac{H_{n}}{\beta}\right)\left(\frac{\kappa_{v}\alpha}{1+\alpha}\right)^{2}\left(\frac{100}{g_{\ast}}\right)^{1/3}v_{w}$
(51)
$\displaystyle\times\left(\frac{f}{f_{sw}}\right)^{3}\left(\frac{7}{4+3(f/f_{\textrm{sw}})^{2}}\right)^{7/2}\Upsilon(\tau_{sw})\
,$

where $v_{w}$ is the wall velocity, and we take $v_{w}=c_{s}=\sqrt{1/3}$ with
$c_{s}$ being the sound velocity. $f_{\text{sw}}$ is the present peak
frequency of the spectrum,

$f_{\textrm{sw}}\ =\ 1.9\times
10^{-5}\frac{1}{v_{w}}\left(\frac{\beta}{H_{n}}\right)\left(\frac{T_{n}}{100\textrm{GeV}}\right)\left(\frac{g_{\ast}}{100}\right)^{1/6}\textrm{Hz}\,.$
(52)

The $\kappa_{v}$ is the fraction of latent heat transformed into the kinetic
energy of the fluid 1004.4187 ,

$\kappa_{v}\simeq\frac{\alpha^{2/5}}{0.017+(0.997+\alpha)^{2/5}}.$ (53)

The suppression factor of Eq. (51) 2007.08537

$\Upsilon(\tau_{sw})=1-\frac{1}{\sqrt{1+2\tau_{sw}H_{n}}},$ (54)

appears due to the finite lifetime $\tau_{sw}$ of the sound waves 2003.07360 ;
2003.08892 ,

$\tau_{sw}=\frac{\tilde{v}_{W}(8\pi)^{1/3}}{\beta\bar{U}_{f}},~{}~{}\bar{U}^{2}_{f}=\frac{3}{4}\frac{\kappa_{v}\alpha}{1+\alpha}.$
(55)

| $\lambda_{SX}$ | $\lambda_{X}$ | $\lambda_{s}$ | $k_{1}k_{2}$ | $s_{L}$ | $s_{R}$ | $g_{Z}\prime$ |  $m_{X_{I}}$(GeV)
---|---|---|---|---|---|---|---|---
BP1 |  0.298 |  1.478 |  0.067 |  0.0061 |  0.113 |  0.254 |  0.481 |  302.9
BP2 |  0.304 |  0.337 |  0.017 |  -0.0164 |  -0.280 |  -0.374 |  0.319 |  319.2

|  $m_{X_{R}}$(GeV) |  $m_{E_{1}}$(GeV) |  $m_{E_{2}}$(GeV) |  $m_{Z}\prime$(GeV) |  $m_{S}$(GeV) |  $T_{n}$(GeV) |  $\xi_{n}$
---|---|---|---|---|---|---|---
BP1 |  410.2 |  357.4 |  454.6 |  308.4 |  117.1 |  214.5 |  1.120
BP2 |  331.4 |  363.8 |  469.0 |  312.3 |  91.6 |  217.4 |  1.571

Table 2: Input and output parameters for the BP1 and BP2. Here
$\xi_{n}=\frac{<S>}{T_{n}}$ denotes the strength of the
$U(1)_{L_{\mu}-L_{\tau}}$ breaking FOPT at the nucleation temperature $T_{n}$.
Figure 4: Gravitational wave spectra for the BP1 and BP2.

We calculate GW spectra for thousands of parameter points accommodating the
muon $g-2$, DM relic density, and the exclusion limits of the LHC direct
searches, and find that all the peak strengths are below the sensitivity curve
of BBO. We pick out two benchmark points (BPs) whose the peak strengths of GW
spectra exceed the sensitivity curve of U-DECIGO. The corresponding input and
output parameters for the two BPs are given in Table 2, and the GW spectra are
shown along with expected sensitivities of various future interferometer
experiments in Fig. 4. For the two BPs, the nucleation temperature of the
$U(1)_{L_{\mu}-L_{\tau}}$ breaking phase transition $T_{n}$ is around 200 GeV.
The lowest peak frequency is 0.003 HZ from BP1, and the highest peak frequency
is 0.2 HZ from BP2.

##  VI Conclusion

We study an extra $U(1)_{L_{\mu}-L_{\tau}}$ gauge symmetry extension of the
standard model by considering the dark matter, muon $g-2$ anomaly, the
$U(1)_{L_{\mu}-L_{\tau}}$ breaking phase transition, GW spectra, and the bound
from the direct detection at the LHC. We obtained the following observations:
(i) A joint explanation of the dark matter relic density and muon $g-2$
anomaly rules out the region where both
$min(m_{E_{1}},m_{E_{2}},m_{N},m_{X_{R}})$ and $min(m_{Z^{\prime}},m_{S})$ are
much larger than $m_{X_{I}}$. (ii) A first order $U(1)_{L_{\mu}-L_{\tau}}$
breaking phase transition can be achieved in the parameter space explaining
the DM relic density and muon $g-2$ anomaly simultaneously, and the
corresponding gravitational wave spectra can reach the sensitivity of
U-DECIGO. (iii) The mass spectra of the vector-like leptons and dark matter
are stringently restricted by the direct searches at the LHC.

## Acknowledgment

This work was supported by the National Natural Science Foundation of China
under grant 11975013.

## References

  * (1) X. G. He, G. C. Joshi, H. Lew and R. R. Volkas, “New-${Z}^{{}^{\prime}}$ phenomenology”, Phys. Rev. D 43, (1991) 22–24.
  * (2) E. Ma, D. P. Roy and S. Roy, Gauged $L_{\mu}-L_{\tau}$ with large muon anomalous magnetic moment and the bimaximal mixing of neutrinos, Phys. Lett. B 525, (2002) 101-106.
  * (3) S. Baek, N. G. Deshpande, X. G. He and P. Ko, Muon anomalous g-2 and gauged $L_{\mu}-L_{\tau}$ models, Phys. Rev. D 64, (2001) 055006.
  * (4) J. Heeck and W. Rodejohann, Gauged $L_{\mu}-L_{\tau}$ Symmetry at the Electroweak Scale, Phys. Rev. D 84, (2011) 075007.
  * (5) K. Harigaya, T. Igari, M. M. Nojiri, M. Takeuchi and K. Tobe, Muon g-2 and LHC phenomenology in the $L_{\mu}-L_{\tau}$ gauge symmetric model, JHEP 03, (2014) 105.
  * (6) W. Altmannshofer, C.-Y. Chen, P. S. Bhupal Dev and A. Soni, Lepton flavor violating Z’ explanation of the muon anomalous magnetic moment, Phys. Lett. B 762, (2016) 389-398.
  * (7) A. Biswas, S. Choubey and S. Khan, Neutrino Mass, Dark Matter and Anomalous Magnetic Moment of Muon in a $U(1)_{L_{\mu}-L_{\tau}}$ Model, JHEP 09, (2016) 147.
  * (8) A. Biswas, S. Choubey and S. Khan, FIMP and Muon ($g-2$) in a U$(1)_{L_{\mu}-L_{\tau}}$ Model, JHEP 02, (2017) 123.
  * (9) H. Banerjee, P. Byakti and S. Roy, Supersymmetric gauged U(1)${}_{L_{\mu}-L_{\tau}}$ model for neutrinos and the muon $(g-2)$ anomaly, Phys. Rev. D 98, (2018) 075022.
  * (10) S. Baek and P. Ko, Phenomenology of $U(1)_{L_{\mu}-L_{\tau}}$ charged dark matter at PAMELA and colliders, JCAP 0910, (2009) 011.
  * (11) M. Das and S. Mohanty, Leptophilic dark matter in gauged $L_{\mu}-L_{\tau}$ extension of MSSM, Phys. Rev. D 89, (2014) 025004.
  * (12) S. Patra, S. Rao, N. Sahoo and N. Sahu, Gauged $U(1)_{L_{\mu}-L_{\tau}}$ model in light of muon $g-2$ anomaly, neutrino mass and dark matter phenomenology, Nucl. Phys. B 917, (2017) 317-336.
  * (13) A. Biswas, S. Choubey, L. Covi and S. Khan, Explaining the 3.5 keV X-ray Line in a ${L_{\mu}-L_{\tau}}$ Extension of the Inert Doublet Model, JCAP 1802, (2018) 002.
  * (14) P. Foldenauer, Light dark matter in a gauged $U(1)_{L_{\mu}-L_{\tau}}$ model, Phys. Rev. D 99, (2019) 035007.
  * (15) A. Crivellin, G. D’Ambrosio and J. Heeck, Phys. Rev. Lett. 114, (2015) 151801.
  * (16) W. Altmannshofer, S. Gori, S. Profumo and F. S. Queiroz, JHEP 12, (2016) 106.
  * (17) C.-H. Chen and T. Nomura, Phys. Lett. B 777, (2018) 420–427.
  * (18) S. Baek, Phys. Lett. B 781, (2018) 376–382.
  * (19) W. Altmannshofer, S. Gori, M. Pospelov and I. Yavin, Phys. Rev. D 89, (2014) 095033.
  * (20) P. Arnan, L. Hofer, F. Mescia and A. Crivellin, JHEP 04, (2017) 043.
  * (21) S. Singirala, S. Sahoo and R. Mohanta, Exploring dark matter, Phys. Rev. D 99, (2019) 035042.
  * (22) P. T. P. Hutauruk, T. Nomura, H. Okada and Y. Orikasa, Phys. Rev. D 99, (2019) 055041.
  * (23) A. Biswas, A. Shaw, JHEP 05, (2019) 165.
  * (24) Z.-L. Han, R. Ding, S.-J. Lin, B. Zhu, Eur. Phys. Jour. C 79, (2019) 1007.
  * (25) A. S. Joshipura, N. Mahajan, K. M. Patel, JHEP 03, (2020) 001.
  * (26) W. Chao, H. Wang, L. Wang and Y. Zhang, “Dark matter, $Z^{\prime}$, vector-like quark at the LHC and $b\to s\mu\mu$ anomaly,” Chin. Phys. C 45, 083105 (2021).
  * (27) LISA Collaboration, H. Audley et al., Laser Interferometer Space Antenna, arXiv:1702.00786.
  * (28) X. Gong et al., Descope of the ALIA mission, J. Phys. Conf. Ser. 610, 012011 (2015).
  * (29) TianQin Collaboration, J. Luo et al., TianQin: a space-borne gravitational wave detector, Class. Quant. Grav. 33, 035010 (2016).
  * (30) K. Yagi and N. Seto, Detector configuration of DECIGO/BBO and identification of cosmological neutron-star binaries, Phys. Rev. D 83, 044011 (2011).
  * (31) H. Kudoh, A. Taruya, T. Hiramatsu, and Y. Himemoto, Detecting a gravitational-wave background with next-generation space interferometers, Phys. Rev. D 73, 064006 (2006).
  * (32) S. Liu and L. Wang, “Spontaneous CP violation electroweak baryogenesis and gravitational wave through multistep phase transitions,” Phys. Rev. D 107, (2023) 115008.
  * (33) W. Altmannshofer, S. Gori, M. Pospelov and I. Yavin, “Neutrino trident production: a powerful probe of new physics with neutrino beams.”, Phys. Rev. Lett. 113, 091801 (2014).
  * (34) P. A. Zyla et al. [Particle Data Group], Review of Particle Physics, PTEP 2020, 083C01 (2020).
  * (35) M. E. Peskin and T. Takeuchi, “New constraint on a strongly interacting Higgs sector”, Phys. Rev. Lett. 65, 964 (1990).
  * (36) M.-C. Chen, S. Dawson, “One-loop radiative corrections to the $\rho$ parameter in the littlest Higgs model”, Phys. Rev. D 70, (2004) 015003.
  * (37) S. K. Garg, C. S. Kim, “Vector like leptons with extended Higgs sector”, arXiv:1305.4712.
  * (38) R. Dermisek, A. Raval, “Explanation of the muon $g\mathbf{-}2$ anomaly with vectorlike leptons and its implications for Higgs decays”, Phys. Rev. D 88, (2013) 013017.
  * (39) J. Kawamura, S. Raby, A. Trautner, “Complete vectorlike fourth family and new ${\mathrm{U}(1)}^{{}^{\prime}}$ for muon anomalies”, Phys. Rev. D 100, (2019) 055030.
  * (40) F. Jegerlehner, A. Nyffeler, “The Muon g-2”, Phys. Rept. 477, (2009) 1–110.
  * (41) B. Abi et al. [Fermilab Collaboration], “Measurement of the Positive Muon Anomalous Magnetic Moment to 0.46 ppm”, Phys. Rev. Lett. 126, (2021) 141801.
  * (42) Muon g-2 Collaboration, “Precise measurement of the positive muon anomalous magnetic moment”, Phys. Rev. Lett. 86, (2001) 2227; “Final Report of the Muon E821 Anomalous Magnetic Moment Measurement at BNL”, Phys. Rev. D 73, (2006) 072003.
  * (43) J. Mott et al., New results from the Muon g-2 experiment at Fermilab, https://indico.fnal.gov/event/60738/.
  * (44) S. Borsanyi et al., Leading-order hadronic vacuum polarization contribution to the muon magnetic momentfrom lattice QCD, arXiv:2002.12347.
  * (45) CMD-3 Collaboration, F. V. Ignatov et al., Measurement of the $e^{+}e^{-}\to Pi^{+}Pi^{-}$ cross section from threshold to 1.2 GeV with the CMD-3 detector, arXiv:2302.08834.
  * (46) Planck Collaboration, “Planck 2015 results. XXVII. The Second Planck Catalogue of Sunyaev-Zeldovich Sources”, Astron. Astrophys. A 27, 594 (2016).
  * (47) G. Belanger, F. Boudjema, A. Pukhov, A. Semenov, “micrOMEGAs-3: A program for calculating dark matter observables”, Comput. Phys. Commun. 185, 960-985 (2014).
  * (48) S. R. Coleman and E. J. Weinberg, Radiative Corrections as the Origin of Spontaneous Symmetry Breaking, Phys. Rev. D 7, 1888 (1973).
  * (49) L. Dolan and R. Jackiw, Symmetry Behavior at Finite Temperature, Phys. Rev. D 9, 3320 (1974).
  * (50) P. B. Arnold and O. Espinosa, The Effective potential and first order phase transitions:Beyond leading-order, Phys. Rev. D 47, 3546 (1993) [Erratum: Phys. Rev. D 50, 6662 (1994)].
  * (51) R. R. Parwani, Resummation in a hot scalar field theory, Phys. Rev. D 45, 4695 (1992).
  * (52) J. M. Cline, K. Kainulainen and M. Trott, “Electroweak Baryogenesis in Two Higgs Doublet Models and B meson anomalies,” JHEP 1111, 089 (2011).
  * (53) J. R. Espinosa, T. Konstandin and F. Riva, “Strong Electroweak Phase Transitions in the Standard Model with a Singlet,” Nucl. Phys. B 854 (2012), 592-630.
  * (54) ATLAS Collaboration, “Search for electroweak production of charginos and sleptons decaying into final states with two leptons and missing transverse momentum in $\sqrt{s}=13$ TeV $pp$ collisions using the ATLAS detector”, Eur. Phys. Jour. C 80, (2020) 123.
  * (55) E. Conte, B. Fuks, “Confronting new physics theories to LHC data with MADANALYSIS 5”, Int. J. Mod. Phys. A 33, (2018) 1830027.
  * (56) J. Y. Araz, M. Frank and B. Fuks, “Reinterpreting the results of the LHC with MadAnalysis 5: uncertainties and higher-luminosity estimates”, Eur. Phys. Jour. C 80, (2020) 531.
  * (57) J. Y. Araz, B. Fuks and G. Polykratis, “Simplified fast detector simulation in MADANALYSIS 5”, Eur. Phys. Jour. C 81, 329 (2021).
  * (58) J. Alwall et al., “The automated computation of tree-level and next-to-leading order differential cross sections, and their matching to parton shower simulations”, JHEP 1407, (2014) 079.
  * (59) P. Torrielli and S. Frixione, “Matching NLO QCD computations with PYTHIA using MC@NLO”, JHEP 1004, (2010) 110.
  * (60) J. de Favereau et al. [DELPHES 3 Collaboration], “DELPHES 3, A modular framework for fast simulation of a generic collider experiment”, JHEP 1402, (2014) 057.
  * (61) M. Hindmarsh, S. J. Huber, K. Rummukainen, and D. J. Weir, Numerical simulations of acoustically generated gravitational waves at a first order phase transition, Phys. Rev. D 92, 123009 (2015).
  * (62) M. Maziashvili, JCAP 1006, 028 (2010).
  * (63) H.-K. Guo, K. Sinha, D. Vagie and G. White, Phase Transitions in an Expanding Universe: Stochastic Gravitational Waves in Standard and Non-Standard Histories, JCAP 01, 001 (2021).
  * (64) J. Ellis, M. Lewicki and J. M. No, Gravitational waves from first-order cosmological phase transitions: lifetime of the sound wave source, JCAP 07, 050 (2020).
  * (65) X. Wang, F. P. Huang and X. Zhang, Phase transition dynamics and gravitational wave spectra of strong first-order phase transition in supercooled universe, JCAP 05, 045 (2020).


11institutetext: IIIT Hyderabad, Telangana 500032, India 11email:
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
# Nostradamus: Weathering Worth

Alapan Chaudhuri  0000-0001-4836-9430 Zeeshan Ahmed  0000-0002-2719-7355
Ashwin Rao  0000-0002-1545-2611 Shivansh Subramanian  0000-0001-5647-1790
Shreyas Pradhan  0000-0002-3762-3653 Abhishek Mittal  0000-0001-8987-8386

###### Abstract

Nostradamus, inspired by the French astrologer and reputed seer, is a detailed
study exploring relations between environmental factors and changes in the
stock market. In this paper, we analyze associative correlation and causation
between environmental elements (including natural disasters, climate and
weather conditions) and stock prices, using historical stock market data,
historical climate data, and various climate indicators such as carbon dioxide
emissions. We have conducted our study based on the US financial market,
global climate trends, and daily weather records to demonstrate a significant
relationship between climate and stock price fluctuation. Our analysis covers
both short-term and long-term rises and dips in company stock performances.
Lastly, we take four natural disasters as a case study to observe the effect
they have on people’s emotional state and their influence on the stock market.

###### Keywords:

Financial data analysis Probabilistic models Statistical inference Stock
market analysis.

##  1 Introduction

A financial market [1] is an aggregation of buyers and sellers trading stocks
or their derivatives. A popular saying in the context of trading stocks goes
that people decide the worth of the company. Though it aptly describes the
definition, we often see that such is not the case. Due to various (expected
or unexpected) reasons, companies end up deviating from market trends. Our
goal here is to look into such deviations which might be resulted from climate
conditions, weather or natural calamities and establish their statistical
significance so that we can predict, analyse and appropriately react to such
events in the future.

##  2 Dataset

###  2.1 Data Extractions

We have taken historical stock data for multiple companies using Yahoo Finance
API [2] and ZEPL US Stock Market Data [3]. We have also extracted
sustainability data about each company from Yahoo for general usage and
insights. All of the data is publicly available for free and is also stored on
our GitHub. For environmental Data, we have used KNOEMA Environmental Data
Atlas [4] to get environmental data for multiple countries over many years. We
have got $40+$ parameters related to CO2 emissions, fuel consumption, and Use
of nitrogen, amongst others per country per year. To get climate-related data,
we have used NOAA Climate Data [5], which gives us the temperature,
precipitation and snowfall of a given area per day.

###  2.2 Data Cleaning and Enhancements

Once we have obtained the data required as data-frames, we performed feature
selection to reduce the noise in the data. This made our process more
efficient and decreased space limitations. We have performed two different
analyses: yearly and daily. For yearly data analysis, we have grouped the
weather data, climate data and historical OHLC data by year to aptly represent
the columns. For daily data analysis, we computed Moving Averages, Volume
Weighted Average Prices [6], Uncertainty and 50 Day Standard Deviation.

$VWAP=\frac{\sum_{0}^{n-1}Volume_{i}\times\frac{High+Low+Close}{3}}{\sum_{0}^{n-1}Volume}$
(1)

$SD=\sqrt{\frac{\sum_{0}^{n-1}\mid x-\bar{x}\mid}{n}}$ (2)

We have stored all the data required in separate CSV files to decrease the
runtime. We also created a pipeline which enables the user to fetch all the
data for any stock they require easily.

##  3 Environmental Factors

We hypothesised that there are multiple environmental factors which affect the
government and the people directly, but not the companies so much. Companies
change their policies, which change their revenues, based on reactions from
the government and their customers. In a graph, this would be reflected by an
initial change in trend (based on people’s opinion), which later settles a bit
and then changes based on whether or not government changes their policies.

$r_{xy}=\frac{\sum(x_{i}-\bar{x})(y_{i}-\bar{y})}{\sqrt{\sum(x_{i}-\bar{x})^{2}\sum(y_{i}-\sum\bar{y})^{2}}}$
(3)

Recent work in this direction includes studies by Castro et al. 2021 [7] and
Zahid et al. 2022 [8]. They looked majorly into the effects on the European
financial markets and South Asian financial markets respectively.

###  3.1 Method

We define the value of correlation $r_{xy}$ as the ratio of the sum of
deviations divided by the root of the product of the sum of squares of
deviations. A high value of the correlation index (close to 1) indicates that
the environmental factor and the company’s profits share a causal
relationship, whereas a correlation index closer to -1 represents an inverse
of such a relationship.

Figure 1: Represents correlation between various companies’ stock prices vs
CO2 emission in the USA over the years (view full-sized image in our GitHub
repository).

###  3.2 Inferences

Multiple stocks show a different correlation between CO2 emissions and stock
prices, and each of them signifies important information related to their
sector.

####  3.2.1 Strong Correlations

A high positive value of correlation usually means that there is an
interdependence (Strong Correlation) between a stock and an environmental
factor. Consider the stock (British Petroleum Company), which is an oil and
gas company. The company’s production has a direct effect on CO2 emissions.
Hence, we can confidently infer that CO2 emission values are affected by BP’s
stock.

Figure 2: Shows the correlation between open prices of BP over the years vs
CO2 emissions in the USA. Correlation value as high as 0.8159 was observed.

####  3.2.2 Causation versus Correlation

It is not always true that a high value of correlation means that there is an
interdependence between stock prices and CO2 emissions. Consider the stock of
AAPL (Apple) and the CO2 emissions from gaseous fuel consumption. The
correlation value observed was close to 0.93. We know that Apple (a tech
company) is not dependent on CO2 emissions from gaseous fuel consumption.
Despite that, it has a high correlation with that factor. This is not a result
of dependence between the two things. The high correlation is simply a
coincidence as CO2 emissions from gaseous fuels are rising because of the
rapid population growth and because nuclear and other clean energy sources are
not very prevalent. Hence correlation is not always a result of causation [9].

Figure 3: Shows the correlation between open prices of AAPL over the years vs
CO2 from gaseous fuel consumption in the USA. Correlation value as high as
0.93 was observed.

####  3.2.3 Hidden Correlations

Sometimes a stock and an unrelated environmental factor have a high
correlation. We observed that there are $2$ possible explanations for it:
either it is just a coincidence (as shown above), or it has a hidden
correlation. For example, consider the stock of EOD (Wells Fargo Global) vs
CO2 emissions in the USA. The graph has a high positive correlation value of
0.927. However, there are likely hidden correlations due to dependencies of
carbon emissions on an industry which also determines the prices of companies
which own/invest in the same. This includes large companies such as large
banks and firms. Hence, this high correlation is very likely not a result of
coincidence but a result of purposeful investing/decisions taken by the
company.

Figure 4: Shows the correlation between open prices of EOD over the years vs
CO2 emissions in the USA. Correlation value as high as 0.927 was observed.

####  3.2.4 Low Positive Correlations

A low positive correlation usually means that the stock and the environmental
factor in consideration are independent of each other. Consider the case of
XOM (Exxon Mobil) and Agricultural Methane Emissions, which has a low
correlation of 0.234. Hence, we can most of the time conclude that these two
are independent of each other. In rare cases, though, they could be dependent
and still have a low correlation.

Figure 5: Shows the correlation between open prices of XOM over the years vs
Agricultural Methane Emissions in the USA. Low correlation value of 0.234 was
observed.

####  3.2.5 High Negative Correlations

Highly negative correlation values generally imply an inverse effect between
the company’s production/success and the factor in consideration. Consider the
volume of the stock CVX (Chevron Corporation) and the factor of Terrestrial
and marine protected areas. This has a highly negative correlation of -0.899.
It is fair to assume that when the number of terrestrial and marine protected
areas increases, the volume of the stock CVX, an energy industry, decreases.
Hence, a highly negative correlation value means inverse dependence.

Figure 6: Shows the correlation between open prices of CVX over the years vs
Terrestrial and Marine Protected areas in the USA. High negative correlation
value of -0.899 was observed.

###  3.3 Predictions

Figure 7: US energy-related CO2 emissions

The global emission levels are predicted to be stranded due to growing
awareness about climate change and a noticeable switch to renewable sources as
primary energy sources.

As we earlier saw the high correlation between the carbon emission levels and
the stock price of BP, it can be analyzed that in the coming times, the stock
price of BP and other large oil companies will begin to fall. This also means
that there will be a rise in the stock prices of companies that provide an
alternative source of fuels, such as solar panels and windmills.

###  3.4 Industrial Revolution

The example of correlation and not causation can further be understood by
considering the rise in the population and workforce. The rise in population
meant more consumption of energy, and since we lack any large source of energy
other than fossil fuels, it means that there will be a rise in pollution
levels.

Figure 8: Correlations between markets and carbon emissions.

If the stock history of the industrial revolution times and the carbon
emission levels of the same were available, it could be seen that the rising
workforce and demand directly meant a rise in carbon levels. This is a
correlation due to familiar sources of the rise.

##  4 Natural Disasters

Now, that we have covered aspects of environmental factors, we wanted to
experiment with the effects of sudden calamities on the stock market.
Therefore, we performed case studies on four disasters (both natural and man-
made) occurring in the USA and how prices of certain relevant stocks reacted
to them. There have been extensive previous studies in this direction as well.
Our conclusions in this respect do concur with several of such works including
Tavor & Teitler-Regev 2019 [10], Nguyen & Chaiechi 2021 (who specifically
looked into the Hong Kong Stock Market) [11] and Seetharam 2017 (a relatively
older paper with an extensive statistical report) [12].

###  4.1 California Wildfires

In terms of property damage, 2017 was the most destructive wildfire season on
record in California then, surpassed by only the 2018 season, with 9,560 fires
burning [13]. Throughout 2017, the fires destroyed or damaged more than 10,000
structures in the state (destroyed 9,470, damaged 810).

Figure 9: Stock prices during the 2017-18 California wildfires. The blue bar
represents when the wildfires started.

We can observe an interdependence between the graphs in this case. The stock
of PCG, which is a gas and electric company, was affected right after the
disaster took place. But now consider EIX. Its stock fell sharply sometime
after PCG’s. We know that EIX depends on PCG for its production. This is also
evident through the graph as its stock was not affected directly after the
wildfires but a few days later. This gives clear information about the
dependence of companies on each other.

###  4.2 Texas Storm

In February 2021, the state of Texas suffered a major power crisis [14], which
came about as a result of three severe winter storms sweeping across the
United States. The storms caused a massive electricity generation failure in
the entire state leading to shortages of water, food, and heat. More than 4.5
million homes and businesses were left without power, and at least 210 people
were killed directly or indirectly, with some estimates going as high as 702.
The Texas Grid failure was majorly caused by the inadequately winterized
natural gas equipment. NRG suffered the most, as can be seen in the following
plots. It can also be emphasized that the analysis has to be done on a small
time scale as companies often bounce back from such losses over a more
extended period of time, such as a year.

Figure 10: Stock prices of NRG during the 2021 Texas winter storm. The blue
bar represents when the winter storm was at its peak.

###  4.3 Katrina Hurricane

In August 2005, Katrina [15] was a category 5 Atlantic hurricane that caused
over 1,800 deaths and $125 billion in damage. This damage was mainly focused
on New Orleans and the surrounding areas. Large transportation companies such
as C.H. Robinson faced heavy losses due to the denial of services. Since
hurricanes have a relatively minor effect on companies, they can resume their
services. so the loss is made back.

Figure 11: Stock prices during the 2005 Katrina Hurricane. The blue bar
represents when the hurricane occurred.

As is seen here, companies related to travel and transport profited here due
to the urgent requirement for vehicles for relief and rescue.

###  4.4 9/11

Although 9/11 [16] was not a natural disaster, it is worth considering how it
affected the stock market. The effect is visible on the entire stock market
due to direct or indirect loss. As seen in the bottom graphs, the prices go
down relatively fast after the event. Due to its surprising nature, the
impact’s nature and magnitude were no less than a natural disaster over the
country.

For all these disasters, stocks of companies related to those events were
considered. We observe that AMK, a company related to airlines, suffered huge
losses. Even stock indices suffered gap losses representing the general mood
of the market and traders after the disasters when everything seemed unsure.
In agreement with our hypothesis, the stocks of those companies suffered a
heavy hit after those events.

Figure 12: Stock prices during 9/11. The blue bar represents when the event
occurred.

##  5 Weather

We try to observe a correlation between daily weather and overall daily trades
and obtain a value of $0.0338$, which can be interpreted as statistically
insignificant. The stock data we have used for this analysis is SPY, and the
weather tracked is of New York City. Our result affirms that of Pardo & Valor,
2003 [17] and Wang & Kuang-Hsunshih & Jang, 2017 [18]. We also found no
correlation between trading volume and precipitation, which supports the
conclusion of Wang & C. Lin & J. Lin, 2011 [19]. However, this result does
contradict several much older works such as Saunders, 1993 [20] and
Hirshleifer & Shumway, 2003 [21]. A possible reason behind such disagreement
might lie in the fact that we are dealing with market data from a much more
recent time period. A lot has changed since the 1990s, and early 2000s
especially because electronic trading has made the financial markets much more
efficient, transparent and location-independent [22].

Figure 13: No particular correlation between weather (NYC daily precipitation
and temperature composite in orange) and daily trading (a composite index
representing SPY volatility and volume in blue).

##  6 Conclusion and Future Work

Our analysis categorically characterises correlations between stock prices and
environmental factors. We can conclude with high confidence that even though
short-term changes in daily weather generally have no effects on daily trades
at the present time, relatively unpredictable events, including natural
disasters have a significant toll on stock prices. Interestingly not all
stocks get affected negatively.

Furthermore, we have studied the financial and economic ties associated with
impactful environmental factors and the cyclicity of catastrophes. Events like
hurricanes and heat waves recur periodically. Due to rising climate change,
they will only become more frequent. This means that past analysis of such
events will help us determine the kind of impact companies would be likely to
face. We hope this work would help in being a stepping stone towards more
studies in this direction. A promising way forward would lie in building a
robust predictive model based on analysis which has already been performed.

## References

  * [1] Mele, Antonio: Financial Economics. MIT Press (2022)
  * [2] Application portal interface (API), Yahoo!Finance, https://finance.yahoo.com/quotes/API,Documentation/view/v1/
  * [3] Zepl: US Stock Market Data for Data Science, Snowflake Datasets, https://www.snowflake.com/datasets/zepl-us-stock-market-data-for-data-science/
  * [4] KNOEMA: Environment Data Atlas, environment-related datasets from OECD, UN, World Bank, EDGAR, and other authoritative sources, Snowflake Datasets, https://www.snowflake.com/datasets/knoema-environment-data-atlas/
  * [5] Climate Data Online (CDO), National Centers for Environmental Information, NOAA, https://www.ncei.noaa.gov/cdo-web/
  * [6] Berkowitz, Stephen et al.: The Total Cost of Transactions on the NYSE. The Journal of Finance, 43(1), 97–112 (2012)
  * [7] Castro, Paula et al.: The impact of environmental performance on stock prices in the green and innovative context. Journal of Cleaner Production 320 (2021)
  * [8] Zahid, RM Ammar et al.: Impact of Environmental Fluctuations on Stock Markets: Empirical Evidence from South Asia. Journal of environmental and public health, vol. 2022:7692086 (2022)
  * [9] Altman, N., Krzywinski, M.: Association, correlation and causation. Nature Methods 12, 899–-900 (2015).
  * [10] Tavor T., Teitler-Regev S.: The impact of disasters and terrorism on the stock market. Journal of Disaster Risk Studies, 11(1), 1–8 (2019).
  * [11] Nguyen T., Chaiechi T.: Chapter 2 - The Effects of Natural Disasters on Stock Market Return and Volatility in Hong Kong. Economic Effects of Natural Disasters, 2021, 11–20, Academic Press (2021)
  * [12] Seetharam, I.: Environmental Disasters and Stock Market Performance. Stanford University, https://web.stanford.edu/~ishuwar/Disasters_Stocks_Current.pdf (2017)
  * [13] 2017 California Wildfires, United States Census Bureau, December 2017, https://www.census.gov/topics/preparedness/events/wildfires/2017-ca-wildfires.html
  * [14] 2021 Winter Storm Uri After-Action Review: Findings Report, City of Austin & Travis County, November 4 (2021), https://www.austintexas.gov/sites/default/files/files/HSEM/2021-Winter-Storm-Uri-AAR-Findings-Report.pdf
  * [15] Gibbens S.: Hurricane Katrina, explained. National Geographic (2019) https://www.nationalgeographic.com/environment/article/hurricane-katrina
  * [16] The 9/11 Commission Report, Final Report of the National Commission on Terrorist Attacks Upon the United States, Official Government Edition (2004), https://www.govinfo.gov/content/pkg/GPO-911REPORT/pdf/GPO-911REPORT.pdf
  * [17] Pardo A., Valor E.: ”Spanish Stock Returns: Where is the Weather Effect?”. European Financial Management 9(1), 117–126 (2003)
  * [18] Wang, YH., Kuang-Hsunshih, K., Jang, Je-Wei: Relationship among Weather Effects, Investors’ Moods and Stock Market Risk: An Analysis of Bull and Bear Markets in Taiwan, Japan and Hong Kong, Panoeconomicus 65(2), 239–253 (2017)
  * [19] Wang, YH., Lin, CT., Lin, J.D.: Does weather impact the stock market? Empirical evidence in Taiwan. Quality & Quantity 46, 695–703 (2012)
  * [20] Saunders, E.: Stock Prices and Wall Street Weather. The American Economic Review, 83(5), 1337–1345 (1993).
  * [21] Hirshleifer, D., Shumway, T.: Good Day Sunshine: Stock Returns and the Weather. The Journal of Finance, 58(3), 1009–1032 (2003).
  * [22] The implications of electronic trading in financial markets, Committee on the Global Financial System, Bank for International Settlements, January 2001, https://www.bis.org/publ/cgfs16.pdf.