# Precise Empirical Determination of Metallicity Dependence of Near-infrared
Period-Luminosity Relations for RR Lyrae Variables

Anupam Bhardwaj INAF-Osservatorio Astronomico di Capodimonte, Via Moiariello
16, 80131 Napoli, Italy  Marcella Marconi INAF-Osservatorio Astronomico di
Capodimonte, Via Moiariello 16, 80131 Napoli, Italy  Marina Rejkuba European
Southern Observatory, Karl-Schwarzschild-Straße 2, 85748, Garching, Germany
Richard de Grijs School of Mathematical and Physical Sciences, Macquarie
University, Balaclava Road, Sydney NSW 2109, Australia  Research Centre for
Astronomy, Astrophysics and Astrophotonics, Macquarie University, Balaclava
Road, Sydney NSW 2109, Australia  Harinder P. Singh Department of Physics and
Astrophysics, University of Delhi, Delhi-110007, India  Vittorio F. Braga
INAF-Osservatorio Astronomico di Roma, via Frascati 33, 00040 Monte Porzio
Catone, Italy  Shashi Kanbur Department of Physics, State University of New
York, Oswego, NY 13126, USA  Chow-Choong Ngeow Graduate Institute of
Astronomy, National Central University, 300 Jhongda Road, 32001 Jhongli,
Taiwan  Vincenzo Ripepi  INAF-Osservatorio Astronomico di Capodimonte, Via
Moiariello 16, 80131 Napoli, Italy  Giuseppe Bono  INAF-Osservatorio
Astronomico di Roma, via Frascati 33, 00040 Monte Porzio Catone, Italy
Dipartimento di Fisica, Università di Roma Tor Vergata, via della Ricerca
Scientifica 1, 00133 Roma, Italy  Giulia De Somma INAF-Osservatorio
Astronomico di Capodimonte, Via Moiariello 16, 80131 Napoli, Italy  INFN-Sez.
di Napoli, Compl. Univ. di Monte S. Angelo, Edificio G, Via Cinthia, 80126
Napoli, Italy  Massimo Dall’Ora  INAF-Osservatorio Astronomico di Capodimonte,
Via Moiariello 16, 80131 Napoli, Italy

(Accepted Feb 9, 2023)

###### Abstract

RR Lyrae variables are excellent population II distance indicators thanks to
their well-defined period-luminosity relations (PLRs) at infrared wavelengths.
We present results of near-infrared (NIR) monitoring of Galactic globular
clusters to empirically quantify the metallicity dependence of NIR PLRs for RR
Lyrae variables. Our sample includes homogeneous, accurate, and precise
photometric data for 964 RR Lyrae variables in 11 globular clusters covering a
large metallicity range ($\Delta\textrm{[Fe/H]}\sim 2$ dex). We derive
$JHK_{s}$ band period-luminosity-metallicity (PLZ) and period-Wesenheit-
metallicity (PWZ) relations anchored using 346 Milky Way field RR Lyrae stars
with Gaia parallaxes, and simultaneously solved for independent distances to
globular clusters. We find a significant metallicity dependence of $\sim 0.2$
mag/dex in $JHK_{s}$ band PLZ and PWZ relations for RR Lyrae stars independent
of the adopted metallicity scale. The metallicity coefficients and the zero-
points of the empirical PLZ and PWZ relations are in excellent agreement with
the predictions from the horizontal branch evolution and pulsation models.
Furthermore, RR Lyrae based distances to our sample of globular clusters are
also statistically consistent with other independent measurements in the
literature. Our recommended empirical $JHK_{s}$ band PLZ relations for RR
Lyrae stars with periods of fundamental mode pulsation ($P_{\textrm{f}}$) are:

$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}M_{J}=-0.44~{}(\pm 0.03)-1.83~{}(\pm
0.02)\log(P_{\textrm{f}})+0.20~{}(\pm
0.02)~{}\textrm{[Fe/H]}~{}~{}(\sigma=0.05~{}\textrm{mag})$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}M_{H}=-0.74~{}(\pm 0.02)-2.29~{}(\pm
0.02)\log(P_{\textrm{f}})+0.19~{}(\pm
0.01)~{}\textrm{[Fe/H]}~{}~{}(\sigma=0.05~{}\textrm{mag})$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}M_{K_{s}}=-0.80~{}(\pm 0.02)-2.37~{}(\pm
0.02)\log(P_{\textrm{f}})+0.18~{}(\pm
0.01)~{}\textrm{[Fe/H]}~{}~{}(\sigma=0.05~{}\textrm{mag})$

††journal: ApJL††software: IDL (Landsman, 1993), Astropy (Astropy
Collaboration et al., 2013)††thanks: Marie-Curie Fellow

##  1 Introduction

RR Lyrae variables are old ($>10$ Gyr), low-mass ($\sim 0.5-0.8M_{\odot}$)
stars that are located at the intersection of the horizontal branch and the
classical instability strip in the Hertzsprung-Russell diagram. These
pulsating stars have long been used as distance indicators thanks to the
correlation of their visual magnitudes with metallicity (see the reviews by
Sandage & Tammann, 2006; Catelan & Smith, 2015; Beaton et al., 2018; Bhardwaj,
2020, and references therein). Since the pioneering work by Longmore et al.
(1986) it has been known that RR Lyrae stars follow a well-defined period-
luminosity relation (PLR) in $K_{s}$-band. Near-infrared (NIR) PLRs offer
advantages over visual magnitude-metallicity relation for distance
determinations due to smaller dependence on interstellar extinction, smaller
variability amplitude, and less sensitivity to evolutionary effects (Bono et
al., 2001; Catelan et al., 2004; Marconi et al., 2015).

Stellar evolution and pulsation models of RR Lyrae stars have predicted a
strong effect of metallicity on NIR PLRs ($\sim 0.18$ mag/dex in $JHK_{s}$,
Catelan et al., 2004; Marconi et al., 2015). The empirical studies however
have reported a wide range for the metallicity coefficient of period-
luminosity-metallicity (PLZ) relations ($0.03-0.23$ mag/dex in $K_{s}$,
Sollima et al., 2006, 2008; Bono et al., 2011; Muraveva et al., 2015, 2018;
Braga et al., 2018; Layden et al., 2019; Neeley et al., 2019; Muhie et al.,
2021; Cusano et al., 2021, and references therein). Most of these empirical
investigations were carried out only in $K_{s}$-band, and very few of them
included all three $JHK_{s}$ filters. Neeley et al. (2019) used Milky Way (MW)
field RR Lyrae stars with parallaxes from Gaia data release 2 (Gaia
Collaboration et al., 2016, 2018) to derive PLZ relations in NIR bands. These
authors found a metallicity coefficient which is consistent with theoretical
predictions but with a four times larger scatter. Bhardwaj et al. (2021b)
showed that the lack of homogeneous NIR photometry and spectroscopic
metallicities limits the precision of PLZ relation for field RR Lyrae stars in
the MW despite the improvements in their parallax measurements with the early
Gaia data release 3 (Gaia Collaboration et al., 2021). The latest Gaia
parallaxes of MW field RR Lyrae stars have also been used to calibrate PLZ
relations at mid-infrared wavelengths based on Wide-field Infrared Survey
Explorer data (e.g., Gilligan et al., 2021; Mullen et al., 2023).

RR Lyrae variables in globular clusters (GCs) of different mean metallicities
provide a unique opportunity to quantify the metallicity coefficient of PLZ
relations. A previous such study was carried out by Sollima et al. (2006) in
$K_{s}$-band using literature photometry of 538 RR Lyrae stars in 16 GCs with
$\Delta\textrm{[Fe/H]}=1.3$ dex, which resulted in a weak constraint on the
metallicity coefficient ($0.08\pm 0.11$ mag/dex). This work aims to use
homogeneous $JHK_{s}$ time-series photometry of RR Lyrae variables in GCs from
our ongoing observational program to empirically calibrate their PLZ relations
and extinction-free period-Wesenheit-metallicity (PWZ) relations. The
empirical calibrations of RR Lyrae PLZ and PWZ relations at NIR wavelengths
will be crucial for stellar astrophysics and distance scale studies based on
population II standard candles. The well-established infrared PLZ relations
for RR Lyrae stars will be particularly important for the James Webb Space
Telescope, upcoming Nancy Grace Roman Space Telescope, and ground-based thirty
meter class telescopes, which will operate predominantly at infrared
wavelengths and enable observations of a large number of RR Lyrae stars in
galaxies beyond the Local Group.

##  2 The Data

Figure 1: Histograms of period distributions (top) of 964 GC RR Lyrae
variables and the distribution of mean metallicities (bottom) for their parent
GCs. The peak and the width of the metallicity distribution represent the mean
metallicity and its scatter from Carretta et al. (2009), respectively. Table
1: Globular clusters with $JHK_{s}$ time-series photometric data.

GC |  | $N$ | Ref.a | [Fe/H]b | E(B-V)c | $\mu_{\textrm{GC}}^{d}$ | $\mu_{\textrm{TW}}^{e}$
---|---|---|---|---|---|---|---
|  | RRab | RRc |  | dex | mag | mag | mag
NGC 6441 | - | 21 | 9 | B22 | $-0.44\pm 0.07$ | 0.47 | $15.52\pm 0.03$ | $15.48\pm 0.02$
NGC 6121 | M4 | 30 | 12 | S14 | $-1.18\pm 0.02$ | 0.35 | $11.34\pm 0.02$ | $11.29\pm 0.02$
NGC 5904 | M5 | 67 | 32 | B23 | $-1.33\pm 0.02$ | 0.03 | $14.37\pm 0.02$ | $14.30\pm 0.02$
NGC 6402 | M14 | 47 | 47 | B23 | $-1.39\pm 0.09$ | 0.60 | $14.80\pm 0.06$ | $14.75\pm 0.02$
NGC 5272 | M3 | 171 | 57 | B20 | $-1.50\pm 0.05$ | 0.01 | $15.04\pm 0.02$ | $15.02\pm 0.02$
NGC 6934 | - | 63 | 12 | B23 | $-1.56\pm 0.09$ | 0.10 | $15.98\pm 0.02$ | $15.95\pm 0.02$
NGC 5139 |  $\omega$ Cen | 83 | 97 | B18 | $-1.64\pm 0.09$ | 0.12 | $13.67\pm 0.02$ | $13.67\pm 0.02$
NGC 7089 | M2 | 22 | 15 | B23 | $-1.66\pm 0.07$ | 0.06 | $15.34\pm 0.02$ | $15.31\pm 0.03$
NGC 5024 | M53 | 28 | 31 | B21a | $-2.06\pm 0.09$ | 0.02 | $16.34\pm 0.02$ | $16.35\pm 0.03$
NGC 7078 | M15 | 48 | 61 | B21b | $-2.33\pm 0.02$ | 0.10 | $15.15\pm 0.02$ | $15.14\pm 0.03$
NGC 6341 | M92 | 8 | 3 | D05 | $-2.35\pm 0.05$ | 0.02 | $14.65\pm 0.02$ | $14.65\pm 0.04$

Note. — aReferences: B22 - Bhardwaj et al. (2022), S14 - Stetson et al.
(2014), B20 - Bhardwaj et al. (2020), B18 - Braga et al. (2018), B21a -
Bhardwaj et al. (2021a), B21b - Bhardwaj et al. (2021b), D05 - Del Principe et
al. (2005), B23 - Bhardwaj et al. (in prep.).

bThe mean metallicities were taken from Carretta et al. (2009).

cThe reddening values were adopted from Harris (2010) with an average
uncertainty of 0.02 mag, which represents the scatter in $E(B-V)$ values from
Harris (2010) and Green et al. (2019).

dThe homogeneous distance moduli ($\mu_{\textrm{GC}}$) are from the catalog of
Baumgardt & Vasiliev (2021).

eThe average distance moduli to GCs derived in this work ($\mu_{\textrm{TW}}$)
using NIR PLZ and PWZ relations in the form of equation (1) for the whole
sample of 1310 RR Lyrae stars.

We use NIR photometry of RR Lyrae variables in GCs from our ongoing
observational programs using the WIRCam instrument at the Canada-France-Hawaii
Telescope (CFHT) and the FLAMINGOS-2 instrument at the Gemini South
telescope111 Interested readers are referred to Bhardwaj et al. (2020) and
Bhardwaj et al. (2022) for details regarding the observing strategy, the data
reduction, and analysis of the data from CFHT and Gemini telescopes,
respectively.. All the photometric data for GCs have been calibrated to Two
Micron All Sky Survey (2MASS, Skrutskie et al., 2006) photometric system.
Time-series photometry has been published for NGC 5272, NGC 5024, NGC 7078,
and NGC 6441 (Bhardwaj et al., 2020, 2021b, 2021a, 2022), and will be
presented for NGC 5904, NGC 6402, NGC 6934, and NGC 7089 in a subsequent study
(Bhardwaj et al. in prep.). We also included homogeneous $JHK_{s}$ mean
magnitudes for RR Lyrae variables in NGC 6341 (Del Principe et al., 2005), NGC
6121 (Stetson et al., 2014), and NGC 5139 (Braga et al., 2018). Our sample
consists of 964 stars including 588 fundamental mode (RRab), 345 first-
overtone (RRc), and 31 mixed-mode RR Lyrae (RRd) variables. The top panel of
Fig. 1 shows period distributions of 964 GC RR Lyrae stars and the bottom
panel shows the distribution of mean metallicities for their parent GCs. Table
1 lists the GCs which host RR Lyrae with $JHK_{s}$ time series photometric
data together with their homogeneous distances, mean metallicities, and
reddening values.

Figure 2: Near-infrared PLRs for MW field RR Lyrae variables in $J$ (top), $H$
(middle) and $K_{s}$ (bottom) bands. The solid black and dashed red lines
represent theoretical PLZ relations for [Fe/H]$=-0.4$ dex and $-2.4$ dex,
respectively. The representative error bars on the left display $\pm 3\sigma$
scatter in the theoretical relations (Marconi et al., 2015).

To calibrate the absolute zero-point of NIR PLZ relations, we also included MW
field RR Lyrae stars with Gaia data and $JHK_{s}$ photometry from 2MASS. The
initial sample of $\sim 400$ stars with positions, metallicities, and
reddening values from the literature was taken from Muraveva et al. (2018).
This sample was cross-matched, within a $1\arcsec$ search radius, with the
2MASS catalog and the Gaia distances and RR Lyrae variability catalogs
(Bailer-Jones et al., 2021; Clementini et al., 2022). The combined sample of
382 stars was further restricted to have metallicity range
($-0.3<\textrm{[Fe/H]}<-2.5$ dex) similar to that of cluster variables and
have 2MASS quality flag ‘A’ in $JHK_{s}$ filters (346 stars). We did not apply
random-phase corrections to single-epoch 2MASS magnitudes because no proper
quantification of amplitude ratios or phase lags is available between $G$ and
$JHK_{s}$ bands, and these observations are not contemporaneous. A comparison
of intensity-averaged magnitudes of 342 RR Lyrae stars in GCs with single-
epoch 2MASS magnitudes (with quality flag ‘A’ in $JHK_{s}$) results in a
median offset of $0.026/0.015/0.051$ mag and a scatter of $0.13/0.11/0.12$ mag
in $J/H/K_{s}$. These systematic offsets were propagated to the photometric
uncertainties on 2MASS magnitudes of MW field stars.

Fig. 2 displays $JHK_{s}$ PLRs for MW field RR Lyrae stars using the final
sample of 346 field variables (321 RRab and 25 RRc stars). The 2MASS $JHK_{s}$
magnitudes were corrected for reddening assuming the reddening law of Cardelli
et al. (1989) with an $R_{V}=3.23$, and the total-to-selective absorption
ratios in NIR filters as: $A_{J/H/K_{s}}=0.94/0.58/0.39E(B-V)$. The periods
and classification of RR Lyrae variables from the Gaia catalog (Clementini et
al., 2022), literature reddening estimates (Muraveva et al., 2018) and the
Gaia based distances (Bailer-Jones et al., 2021) were used to obtain absolute
magnitudes for individual stars. A large scatter in the PLZ relations
($0.20/0.17/0.16$ mag in $J/H/K_{s}$) for MW field stars is seen in Fig. 2 as
compared to the theoretical relations ($0.06/0.04/0.04$ mag in $J/H/K_{s}$,
Marconi et al., 2015), which is not only due to random-phase magnitudes but
also due to parallax uncertainties, heterogeneous metallicity and reddening
estimates. Nevertheless, the slopes of NIR PLRs for RR Lyrae variables in GCs
are well constrained (Bhardwaj, 2022), which will allow an accurate
determination of absolute zero-points of PLZ and PWZ relations at $JHK_{s}$
wavelengths using MW field stars. This analysis will use the combined sample
of 1310 RR Lyrae stars from GCs (964 sources) and MW field (346 sources).

##  3 Linear regression formalism

Our goal is to simultaneously solve for a common slope and a metallicity
coefficient, and a separate zero-point for each GC. This approach does not
require a prior assumption of an independent distance to a given GC, which can
be obtained simultaneously given an absolute zero-point based on field RR
Lyrae stars. In the mathematical form, the extinction-corrected apparent
magnitude or the Wesenheit magnitude of the $j^{\textrm{th}}$ RR Lyrae star in
$i^{\textrm{th}}$ GC at wavelength $\lambda$ is defined as:

$\displaystyle
m_{\lambda_{i,j}}=\mu_{i}+a_{\lambda}+b_{\lambda}\log(P_{i,j})+c_{\lambda}\textrm{[Fe/H]}_{i,j},$
(1)

where, $P_{i,j}$ is the period in days and [Fe/H]i,j is the metallicity of RR
Lyrae stars. For a given cluster, the mean metallicity from Carretta et al.
(2009), listed in Table 1, was adopted for all RR Lyrae members. The term
($\mu_{i}+a_{\lambda}$) also gives the apparent zero-point of the PLZ/PWZ
relation in a given cluster. In the case of field variables, the absolute
magnitude of a given RR Lyrae in the above equation is:

$\displaystyle
M_{\lambda_{j}}=a_{\lambda}+b_{\lambda}\log(P_{j})+c_{\lambda}\textrm{[Fe/H]}_{j}.$
(2)

To solve equations (1) and (2) simultaneously, we setup a matrix formalism
where vector $Y$ includes all the magnitude measurements. The set of equations
are given by a matrix $L$ of $n_{Y}$ rows and $n_{Q}$ columns. Here, $n_{Y}$
is the total number of magnitude measurements in $Y$ and $n_{Q}$ is the number
of free parameters (11 GC distances, zero-point, slope, and metallicity
coefficient) in the vector $Q$. The $\chi^{2}$ statistics is given as
$\chi^{2}=(Y-LQ)^{T}C^{-1}(Y-LQ)$, where $C$ is the error matrix. The maximum-
likelihood parameters are obtained using $Q=(L^{T}C^{-1}L)^{-1}L^{T}C^{-1}Y$,
while the standard errors on the free parameters come from the covariance
matrix $(L^{T}C^{-1}L)^{-1}$. This setup of matrix equations is similar to
those described in detail in Riess et al. (2016, in the appendix). This matrix
formalism is employed to obtain free parameters $Q$, which are used to compute
residuals of the best-fitting linear-regression. This process is iterated by
removing the single largest outlier in each iteration until all residuals are
within $3\sigma$, where $\sigma$ is the standard deviation of the global fit.
We created $10^{3}$ random realizations by varying the magnitudes, adopted
mean-metallicity and reddening values within their uncertainties to obtain
final coefficients and their associated uncertainties.

##  4 Period-Luminosity-Metallicity relations

Figure 3: Near-infrared PLZ and PWZ relations for RR Lyrae variables in GCs in
$JHK_{s}$ bands. In each case, the color-coded solid lines in the top panels
correspond to [Fe/H]=[$-0.5,-1.5,-2.5$] dex. The solid line in the bottom
panels shows the best-fitting linear regression while the dashed line displays
the predicted relation from Marconi et al. (2015). The average residuals
($\Delta J,~{}\Delta H,~{}\Delta K_{s}$) in a period bin of $\log(P)=0.02$
days are shown for visualization purposes. The data points are color-coded
according to the color bar.

We used $JHK_{s}$ mean magnitudes for 964 RR Lyrae variables in GCs and 2MASS
single epoch magnitudes for 346 MW field stars jointly to derive the PLZ and
PWZ relations. Three different samples of RR Lyrae stars were considered
separately: RRab, RRc+RRd, and the combined sample of all 1310 RR Lyrae
variables. The mixed-mode RRd variables have the dominant first-overtone
periods as seen in Fig. 1. Therefore, the first-overtone periods of RRc and
RRd stars were fundamentalized using the equation:
$\log(P_{f})=\log(P_{o})+0.127$ (e.g., Coppola et al., 2015; Braga et al.,
2022) for the combined sample of all 1310 RR Lyrae stars.

The intensity-averaged $JHK_{s}$ magnitudes for variables in GCs were
corrected for reddening assuming the previously mentioned reddening law of
Cardelli et al. (1989) and absorption ratios, and using $E(B-V)$ values
provided by Harris (2010) listed in Table 1. The extinction-free Wesenheit
magnitudes (Madore, 1982) in NIR filters were defined as in Bhardwaj et al.
(2016):
$W_{J,H}=H-1.83(J-H),~{}W_{J,K_{s}}=K-0.69(J-K_{s}),~{}W_{H,K_{s}}=K-1.92(H-K_{s})$.

We simultaneously solved for individual GC distances and derived PLZ relations
for the global sample of 1310 RR Lyrae stars. In the case of RRc variables, we
also added 52 RRc models with $0.25<P<0.45$ days (Marconi et al., 2015) to the
calibrator sample because there are only 25 MW field RRc stars. These RRc
models allowed us to break the degeneracy between the metallicity coefficient
and the zero-points. Table 2 lists the results of the best-fitting linear
regression formalism discussed in Section 3. Fig. 3 displays the PLZ/PWZ
relations for the combined sample of RR Lyrae variables in GCs. The
metallicity effect on RR Lyrae magnitudes is distinctly evident in the top
panels showing that the metal poor stars are brighter. Tighter PLZ/PWZ
relations are seen in the bottom panels where the metallicity term is taken
into account, and the residuals do not exhibit any significant trend as a
function of period or metallicity. These empirical metallicity coefficients
are the most precisely determined in $JHK_{s}$ filters, and are in excellent
agreement with the horizontal branch evolution and pulsation model
predictions. There is a hint of increasing metallicity dependence at shorter
wavelengths, but the differences are not statistically significant. We did not
find any statistically significant changes in the PLZ relations for different
metallicity scales and reddening values as discussed in the appendix A.

Table 2: NIR period-luminosity-metallicity and period-Wesenheit-metallicity
relations.

Band | Type | a | b | c | $\sigma$ | $N_{\textrm{f}}$ | $\chi^{2}_{\textrm{dof}}$ | $\Delta\mu$ | $\Delta\textrm{\bf a}$ | $\Delta\textrm{\bf b}$ | $\Delta\textrm{\bf c}$
---|---|---|---|---|---|---|---|---|---|---|---
$J$ | RRL | $-0.44\pm 0.03$ | $-1.83\pm 0.02$ | $0.20\pm 0.02$ | 0.05 | 1049 | 1.03 | $0.02\pm 0.05$ | $-0.06$ | $-0.07$ | $-0.02$
$H$ | RRL | $-0.74\pm 0.02$ | $-2.29\pm 0.02$ | $0.19\pm 0.01$ | 0.05 | 1037 | 0.98 | $-0.01\pm 0.04$ | $-0.02$ | $0.06$ | $-0.01$
$K_{s}$ | RRL | $-0.80\pm 0.02$ | $-2.37\pm 0.02$ | $0.18\pm 0.01$ | 0.05 | 1077 | 1.02 | $-0.01\pm 0.03$ | $-0.02$ | $0.12$ | $-0.00$
$W_{J,H}$ | RRL | $-1.20\pm 0.02$ | $-2.96\pm 0.02$ | $0.18\pm 0.01$ | 0.07 | 989 | 1.07 | $-0.05\pm 0.04$ | $0.01$ | $0.20$ | $0.02$
$W_{J,K_{s}}$ | RRL | $-1.06\pm 0.03$ | $-2.73\pm 0.02$ | $0.16\pm 0.02$ | 0.06 | 1096 | 0.99 | $-0.04\pm 0.03$ | $0.01$ | $0.23$ | $0.02$
$W_{H,K_{s}}$ | RRL | $-0.93\pm 0.02$ | $-2.45\pm 0.02$ | $0.16\pm 0.01$ | 0.07 | 1033 | 1.17 | $-0.01\pm 0.04$ | $-0.01$ | $0.13$ | $0.02$
$J$ | RRab | $-0.46\pm 0.02$ | $-1.97\pm 0.04$ | $0.20\pm 0.01$ | 0.05 | 651 | 0.83 | $0.03\pm 0.05$ | $-0.05$ | $-0.01$ | $0.01$
$H$ | RRab | $-0.74\pm 0.02$ | $-2.28\pm 0.03$ | $0.19\pm 0.01$ | 0.05 | 667 | 0.88 | $-0.01\pm 0.04$ | $-0.02$ | $0.04$ | $-0.00$
$K_{s}$ | RRab | $-0.81\pm 0.03$ | $-2.39\pm 0.03$ | $0.18\pm 0.02$ | 0.05 | 709 | 0.99 | $-0.02\pm 0.04$ | $-0.01$ | $0.12$ | $0.00$
$W_{J,H}$ | RRab | $-1.19\pm 0.03$ | $-2.91\pm 0.04$ | $0.18\pm 0.02$ | 0.07 | 632 | 0.88 | $-0.06\pm 0.04$ | $0.02$ | $0.24$ | $0.02$
$W_{J,K_{s}}$ | RRab | $-1.07\pm 0.03$ | $-2.69\pm 0.04$ | $0.16\pm 0.02$ | 0.06 | 714 | 0.84 | $-0.05\pm 0.05$ | $0.04$ | $0.20$ | $0.03$
$W_{H,K_{s}}$ | RRab | $-0.93\pm 0.03$ | $-2.41\pm 0.04$ | $0.16\pm 0.01$ | 0.07 | 680 | 0.99 | $0.00\pm 0.06$ | $0.01$ | $0.07$ | $0.01$
$J$ | RRc+RRd | $-0.89\pm 0.04$ | $-2.14\pm 0.05$ | $0.17\pm 0.03$ | 0.05 | 437 | 1.02 | $0.03\pm 0.04$ | $-0.18$ | $-0.32$ | $-0.02$
$H$ | RRc+RRd | $-1.23\pm 0.04$ | $-2.51\pm 0.05$ | $0.17\pm 0.03$ | 0.05 | 428 | 0.93 | $0.04\pm 0.02$ | $-0.08$ | $-0.19$ | $-0.01$
$K_{s}$ | RRc+RRd | $-1.30\pm 0.04$ | $-2.60\pm 0.05$ | $0.17\pm 0.03$ | 0.05 | 435 | 1.03 | $0.04\pm 0.03$ | $-0.08$ | $-0.12$ | $-0.02$
$W_{J,H}$ | RRc+RRd | $-1.72\pm 0.04$ | $-2.92\pm 0.06$ | $0.18\pm 0.03$ | 0.06 | 392 | 1.08 | $0.09\pm 0.05$ | $0.00$ | $-0.16$ | $-0.01$
$W_{J,K_{s}}$ | RRc+RRd | $-1.55\pm 0.04$ | $-2.88\pm 0.05$ | $0.18\pm 0.03$ | 0.05 | 422 | 1.07 | $0.06\pm 0.03$ | $-0.03$ | $-0.02$ | $-0.02$
$W_{H,K_{s}}$ | RRc+RRd | $-1.37\pm 0.04$ | $-2.71\pm 0.05$ | $0.18\pm 0.03$ | 0.06 | 410 | 1.22 | $0.07\pm 0.04$ | $-0.11$ | $-0.04$ | $-0.03$

Note. — The zero-point ($a$), slope ($b$), metallicity coefficient ($c$),
dispersion ($\sigma$), number of stars ($N_{f}$), and the chi-square per
degree of freedom ($\chi^{2}_{\textrm{dof}}$) of the final PLR fits are
tabulated. $\Delta\mu$ represents the mean and the standard deviation of
residuals between distance moduli from this work and those from Baumgardt &
Vasiliev (2021). The last three columns show the difference in the zero-
points, the slopes, and the metallicity coefficients between predicted
(Marconi et al., 2015) and empirical relations derived in this study. In the
case of RRc+RRd stars only, the calibrator sample of MW field stars also
includes 52 models from Marconi et al. (2015) as mentioned in the text.

Figure 3 also shows predicted PLZ/PWZ relations from Marconi et al. (2015,
corresponding to [Fe/H]$~{}=-0.5,-1.5$ and $-2.5$ dex) which suggest a good
agreement with our empirical relations. The last three columns of Table 2 list
the differences in the slopes, zero-points, and the metallicity coefficients
between theoretical (Marconi et al., 2015) and empirical PLZ/PWZ relations
derived in this work. While the metallicity coefficients are very similar,
there are some differences in the slopes and zero-points for some specific
sample and/or filter combinations. In particular, the slope and zero-point of
RRc stars in $J$-band are significantly shallower and fainter than the model
predictions, respectively. We note that the error in the predicted slope for
RRc stars in $J$-band is $0.08$ mag (Marconi et al., 2015), and therefore, the
difference of $-0.32$ mag in the slopes between predicted and empirical
relations is just beyond $3\sigma$ of their combined uncertainties. Similarly,
the largest difference in zero-point ($-0.18\pm 0.06$ mag) also occurs for RRc
stars in the $J$-band, but it is within $\pm 0.05$ mag for most sample and/or
filter combinations.

###  4.1 Distance moduli to globular clusters

Table 1 also lists the average distance moduli to the GCs derived after
solving equations (1) and (2), for PLZ and PWZ relations for the sample of all
1310 RR Lyrae stars. Fig. 4 in the appendix B shows no obvious trend in the
difference between the distance moduli derived in this work and those from
Baumgardt & Vasiliev (2021) as a function of metallicity. While the residual
offset in average distance moduli is statistically small for all RR Lyrae
sample ($-0.026\pm 0.025$ mag), the RR Lyrae distances to GCs are
systematically smaller using empirical calibrations for all/RRab sample and
larger using semi-empirical calibrations for RRc+RRd sample (see $\Delta\mu$
in Table 2). Bhardwaj et al. (2021b) also found a similar result with absolute
zero-points of $-0.76/-0.84$ mag in $K_{s}$ using empirical/theoretical
calibrations. The empirical calibrations use Gaia distances taking into
account the parallax zero-point offset from Lindegren et al. (2021). A fainter
absolute zero-point of empirical PLZ/PWZ relation, which results in smaller
distance moduli, suggests an over-correction in parallaxes because a smaller
parallax zero-point offset than the Lindegren et al. (2021) correction will
result in a brighter zero-point for RR Lyrae stars (see Bhardwaj et al.,
2021b). A large scatter in the PLZ relations for MW field stars (Fig. 2)
precludes us from quantifying this small offset in parallaxes in the latest
Gaia data. Comparison of distance moduli listed in Table 1 also suggests that
the difference is larger for metal-rich GCs while the metal-poor GCs show
excellent agreement. This difference can be attributed to a small number of
both RR Lyrae stars and GCs for [Fe/H] $<-1.2$ dex (see Fig. 1), which could
affect the accuracy of the apparent zero-points of PLZ relations, and in turn
distance moduli, for the relatively metal-rich GCs.

###  4.2 Validation of empirical period-luminosity-metallicity relations

To validate the accuracy of our empirically derived PLZ and PWZ relations, we
determined a distance to the LMC using available NIR photometry in the
literature. We use $\sim 18500$ RR Lyrae stars (13611 RRab, 4848 RRc) from
Cusano et al. (2021) with $JK_{s}$ photometry from the VMC survey together
with their photometric metallicity and reddening estimates. Our empirical
calibrations result in an average distance modulus to the LMC of $18.54\pm
0.13$, $18.57\pm 0.12$, and $18.60\pm 0.15$ mag in $K_{s}$ band using the
whole sample, RRab, and RRc+RRd stars, respectively. The quoted values and
their uncertainties represent the peak and the width of the Gaussian
distribution taking into account the uncertainties on the metallicity scale
and on individual metallicities. When using theoretical relations of Marconi
et al. (2015), we obtained a distance of $18.60\pm 0.13$, and $18.59\pm 0.15$
mag for all/RRab and RRc+RRd samples, respectively, which reflects the
difference in the slopes and zero-points between empirical and theoretical
calibrations. The difference in the LMC distance modulus with the percent-
level precise value of 18.477$\pm$0.026 mag (Pietrzyński et al., 2019) is
likely due to the differences in the calibrator and target PLZ relations, and
due to uncertainties in photometric metallicities and reddening, and
geometrical effects, which contribute to three times larger scatter in the PLZ
relations for RR Lyrae stars in the LMC.

##  5 Conclusions

We presented the results of a NIR monitoring program of GCs to quantify the
metallicity dependence on RR Lyrae PLZ and PWZ relations. The largest sample
of 964 cluster RR Lyrae stars with homogeneous $JHK_{s}$ mean magnitudes based
on time-series data in 11 GCs of different mean metallicities was used in this
analysis. When anchored using additional 346 MW field RR Lyrae stars with Gaia
parallaxes, we obtained purely empirical PLZ and PWZ relations that agree well
with the theoretical predictions based on horizontal branch evolution and
pulsation models. The metallicity coefficient of NIR PLZ and PWZ relations is
$\sim 0.2$ mag/dex, in excellent agreement with the predicted relations. We
simultaneously obtained distance moduli to individual GCs and found those to
be consistent with most recent and independent measurements in the literature
based on Gaia data. While the metallicity coefficients and the slopes are
firmly constrained, the absolute zero-points of PLZ/PWZ relations are
sensitive, at the percent-level precision, to the uncertainties in parallaxes,
parallax zero-point offset, lack of homogeneous time-series photometry and
spectroscopic metallicities for MW field stars. Furthermore, any possible
variations in the slopes of the underlying PLZ/PWZ relations in different
stellar environments may also affect the zero-point, and this basic assumption
of universality of the PLRs will be investigated in detail in a future study.
Nevertheless, our accurate and precise empirical calibrations of PLZ/PWZ
relations will be crucial for RR Lyrae based distance measurements and the
calibration of the distance scale based on population II standard candles in
the era of space telescopes and next-generation ground-based facilities
operating mainly at infrared wavelengths.

We thank the anonymous referee for the constructive referee report that helped
improve the manuscript. This project has received funding from the European
Union’s Horizon 2020 research and innovation programme under the Marie
Skłodowska-Curie grant agreement No. 886298. This research was supported in
part by the Australian Research Council Centre of Excellence for All Sky
Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013.
HPS acknowledges grant 03(1428)/18/EMR-II from Council of Scientific and
Industrial Research (CSIR), India. CCN thanks the funding from the National
Science and Technology Council (Taiwan) under the contract
109-2112-M-008-014-MY3. Access to the CFHT was made possible by the Institute
of Astronomy and Astrophysics, Academia Sinica. CFHT-WIRCam, Gemini
FLAMINGOS-2

##  Appendix A Impact of different metallicity scales and reddening values

We tested for possible variations in the PLZ relations due to different
metallicity scales and reddening values in the literature using the full
sample of 1310 RR Lyrae stars. Our calibrator sample of 346 RR Lyrae stars was
adopted from Muraveva et al. (2018) which was based on the homogenized
compilation of spectroscopic and photometric metallicities by Dambis et al.
(2013) in Zinn & West (1984) metallicity scale. Crestani et al. (2021)
published high resolution spectroscopic metallicities for 208 MW field RR
Lyrae stars. There are 166 RR Lyrae stars in the Crestani et al. (2021) sample
after restricting that to the same metallicity range and 2MASS quality flag
criteria as in Section 2. The coefficients of the PLZ relations derived using
these 166 stars together with 964 GCs variables are statistically consistent
with those listed in Table 2 within $1\sigma$ of their quoted uncertainties.
This was expected considering a small median offset of
$\Delta$[Fe/H]$~{}=0.09\pm 0.02$ dex for 133 common stars in the two
calibrator samples of MW field stars.

In the case of mean metallicities of GCs, the two sets of independent
estimates were obtained from the catalogs of Harris (2010) and Dias et al.
(2016). The mean difference in the [Fe/H] values of 11 GCs is $0.03$ dex
between the estimates from Harris (2010) and Carretta et al. (2009). This
difference in the [Fe/H] values is $0.08$ dex between the measurements of Dias
et al. (2016) and Carretta et al. (2009). Given these small mean differences
with our adopted metallicities, the coefficients of PLZ relations were found
to be in agreement with the values listed in Table 2 within $0.5\sigma$ of
their uncertainties, when using the aforementioned two independent sets of
mean [Fe/H] values.

To investigate the impact of the adopted extinction law and reddening values,
we obtained $E(B-V)$ values from the Baystar19 3D dust map (Green et al.,
2019). The Baystar19 map does not cover NGC 6441 or NGC 5139, for which
$E(B-V)$ values were adopted from Schlegel et al. (1998). The $E(B-V)$ values
from these maps are in the same units and the corresponding absorption ratios
in $JHK_{s}$ were taken from Green et al. (2019, Table 1). After applying
these independent reddening corrections, we derive $JHK_{s}$ band PLZ
relations which exhibit the same slope as listed in Table 2. The metallicity
coefficients increase by 0.01 dex in all three NIR filters and the zero-points
are fainter by 0.04/0.02/0.02 mag in $J/H/K_{s}$ bands. These differences in
the $J$ and $HK_{s}$ band zero-points are within $1.4\sigma$ and $1\sigma$ of
their associated uncertainties, respectively.

##  Appendix B Difference in the distance moduli to globular clusters

Fig. 4 displays the difference in the distance moduli listed in the last two
columns of Table 1.

Figure 4: The difference between the distance moduli ($\mu_{\textrm{TW}}$)
resulting from equation (1) and those from Baumgardt & Vasiliev (2021,
$\mu_{\textrm{GC}}$), as listed in Table 1. The mean and the standard
deviation of residuals are quoted at the bottom left. The dashed and dotted
lines represent the mean and three times the standard deviation of the
residuals, respectively.

## References

  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Bailer-Jones et al. (2021) Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Demleitner, M., & Andrae, R. 2021, AJ, 161, 147, doi: 10.3847/1538-3881/abd806
  * Baumgardt & Vasiliev (2021) Baumgardt, H., & Vasiliev, E. 2021, MNRAS, 505, 5957, doi: 10.1093/mnras/stab1474
  * Beaton et al. (2018) Beaton, R. L., Bono, G., Braga, V. F., et al. 2018, Space Science Reviews, 214, 113, doi: 10.1007/s11214-018-0542-1
  * Bhardwaj (2020) Bhardwaj, A. 2020, Journal of Astrophysics and Astronomy, 41, 23, doi: 10.1007/s12036-020-09640-z
  * Bhardwaj (2022) —. 2022, Universe, 8, 122, doi: 10.3390/universe8020122
  * Bhardwaj et al. (2016) Bhardwaj, A., Kanbur, S. M., Macri, L. M., et al. 2016, AJ, 151, 88, doi: 10.3847/0004-6256/151/4/88
  * Bhardwaj et al. (2022) Bhardwaj, A., Kanbur, S. M., Rejkuba, M., et al. 2022, A&A, 668, A59, doi: 10.1051/0004-6361/202244728
  * Bhardwaj et al. (2020) Bhardwaj, A., Rejkuba, M., de Grijs, R., et al. 2020, AJ, 160, 220, doi: 10.3847/1538-3881/abb3f9
  * Bhardwaj et al. (2021a) Bhardwaj, A., Rejkuba, M., Sloan, G. C., Marconi, M., & Yang, S.-C. 2021a, ApJ, 922, 20, doi: 10.3847/1538-4357/ac214d
  * Bhardwaj et al. (2021b) Bhardwaj, A., Rejkuba, M., de Grijs, R., et al. 2021b, ApJ, 909, 200, doi: 10.3847/1538-4357/abdf48
  * Bono et al. (2011) Bono, G., Dall’Ora, M., Caputo, F., et al. 2011, in RR Lyrae Stars, Metal-Poor Stars, and the Galaxy, ed. A. McWilliam, Vol. 5, 1.  https://arxiv.org/abs/1108.5372
  * Bono et al. (2001) Bono, G., Gieren, W. P., Marconi, M., Fouqué, P., & Caputo, F. 2001, ApJ, 563, 319, doi: 10.1086/323834
  * Braga et al. (2018) Braga, V. F., Stetson, P. B., Bono, G., et al. 2018, AJ, 155, 137, doi: 10.3847/1538-3881/aaadab
  * Braga et al. (2022) Braga, V. F., Fiorentino, G., Bono, G., et al. 2022, MNRAS, 517, 5368, doi: 10.1093/mnras/stac2813
  * Cardelli et al. (1989) Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245, doi: 10.1086/167900
  * Carretta et al. (2009) Carretta, E., Bragaglia, A., Gratton, R., D’Orazi, V., & Lucatello, S. 2009, A&A, 508, 695, doi: 10.1051/0004-6361/200913003
  * Catelan et al. (2004) Catelan, M., Pritzl, B. J., & Smith, H. A. 2004, ApJS, 154, 633, doi: 10.1086/422916
  * Catelan & Smith (2015) Catelan, M., & Smith, H. A. 2015, Pulsating Stars (Wiley-VCH)
  * Clementini et al. (2022) Clementini, G., Ripepi, V., Garofalo, A., et al. 2022, A&A, arXiv:2206.06278, doi: 10.1051/0004-6361/202243964
  * Coppola et al. (2015) Coppola, G., Marconi, M., Stetson, P. B., et al. 2015, ApJ, 814, 71, doi: 10.1088/0004-637X/814/1/71
  * Crestani et al. (2021) Crestani, J., Braga, V. F., Fabrizio, M., et al. 2021, ApJ, 914, 10, doi: 10.3847/1538-4357/abfa2310.48550/arXiv.2104.08113
  * Cusano et al. (2021) Cusano, F., Moretti, M. I., Clementini, G., et al. 2021, MNRAS, 504, 1, doi: 10.1093/mnras/stab901
  * Dambis et al. (2013) Dambis, A. K., Berdnikov, L. N., Kniazev, A. Y., et al. 2013, MNRAS, 435, 3206, doi: 10.1093/mnras/stt1514
  * Del Principe et al. (2005) Del Principe, M., Piersimoni, A. M., Bono, G., et al. 2005, AJ, 129, 2714, doi: 10.1086/430148
  * Dias et al. (2016) Dias, B., Barbuy, B., Saviane, I., et al. 2016, A&A, 590, A9, doi: 10.1051/0004-6361/201526765
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1, doi: 10.1051/0004-6361/201629272
  * Gaia Collaboration et al. (2018) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1, doi: 10.1051/0004-6361/201833051
  * Gaia Collaboration et al. (2021) —. 2021, A&A, 649, A1, doi: 10.1051/0004-6361/202039657
  * Gilligan et al. (2021) Gilligan, C. K., Chaboyer, B., Marengo, M., et al. 2021, MNRAS, 503, 4719, doi: 10.1093/mnras/stab857
  * Green et al. (2019) Green, G. M., Schlafly, E., Zucker, C., Speagle, J. S., & Finkbeiner, D. 2019, ApJ, 887, 93, doi: 10.3847/1538-4357/ab5362
  * Harris (2010) Harris, W. E. 2010, arXiv:1012.3224, arXiv:1012.3224.  https://arxiv.org/abs/1012.3224
  * Landsman (1993) Landsman, W. B. 1993, in ASP Conference Series, Vol. 52, Astronomical Data Analysis Software and Systems II, ed. R. J. Hanisch, R. J. V. Brissenden, & J. Barnes, 246
  * Layden et al. (2019) Layden, A. C., Tiede, G. P., Chaboyer, B., Bunner, C., & Smitka, M. T. 2019, AJ, 158, 105, doi: 10.3847/1538-3881/ab2e10
  * Lindegren et al. (2021) Lindegren, L., Klioner, S. A., Hernández, J., et al. 2021, A&A, 649, A2, doi: 10.1051/0004-6361/202039709
  * Longmore et al. (1986) Longmore, A. J., Fernley, J. A., & Jameson, R. F. 1986, MNRAS, 220, 279, doi: 10.1093/mnras/220.2.279
  * Madore (1982) Madore, B. F. 1982, ApJ, 253, 575, doi: 10.1086/159659
  * Marconi et al. (2015) Marconi, M., Coppola, G., Bono, G., et al. 2015, ApJ, 808, 50, doi: 10.1088/0004-637X/808/1/50
  * Muhie et al. (2021) Muhie, T. D., Dambis, A. K., Berdnikov, L. N., Kniazev, A. Y., & Grebel, E. K. 2021, MNRAS, 502, 4074, doi: 10.1093/mnras/stab219
  * Mullen et al. (2023) Mullen, J. P., Marengo, M., Martínez-Vázquez, C. E., et al. 2023, arXiv e-prints, arXiv:2301.03777.  https://arxiv.org/abs/2301.03777
  * Muraveva et al. (2018) Muraveva, T., Delgado, H. E., Clementini, G., Sarro, L. M., & Garofalo, A. 2018, MNRAS, 481, 1195, doi: 10.1093/mnras/sty2241
  * Muraveva et al. (2015) Muraveva, T., Palmer, M., Clementini, G., et al. 2015, ApJ, 807, 127, doi: 10.1088/0004-637X/807/2/127
  * Neeley et al. (2019) Neeley, J. R., Marengo, M., Freedman, W. L., et al. 2019, MNRAS, 490, 4254, doi: 10.1093/mnras/stz2814
  * Pietrzyński et al. (2019) Pietrzyński, G., Graczyk, D., Gallenne, A., et al. 2019, Nature, 567, 200, doi: 10.1038/s41586-019-0999-4
  * Riess et al. (2016) Riess, A. G., Macri, L. M., Hoffmann, S. L., et al. 2016, ApJ, 826, 56, doi: 10.3847/0004-637X/826/1/56
  * Sandage & Tammann (2006) Sandage, A., & Tammann, G. A. 2006, ARA&A, 44, 93, doi: 10.1146/annurev.astro.43.072103.150612
  * Schlegel et al. (1998) Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525, doi: 10.1086/305772
  * Skrutskie et al. (2006) Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163, doi: 10.1086/498708
  * Sollima et al. (2008) Sollima, A., Cacciari, C., Arkharov, A. A. H., et al. 2008, MNRAS, 384, 1583, doi: 10.1111/j.1365-2966.2007.12804.x
  * Sollima et al. (2006) Sollima, A., Cacciari, C., & Valenti, E. 2006, MNRAS, 372, 1675, doi: 10.1111/j.1365-2966.2006.10962.x
  * Stetson et al. (2014) Stetson, P. B., Braga, V. F., Dall’Ora, M., et al. 2014, PASP, 126, 521, doi: 10.1086/677195
  * Zinn & West (1984) Zinn, R., & West, M. J. 1984, ApJs, 55, 45, doi: 10.1086/190947


# Correlating Gravitational Waves with $W$-boson Mass, FIMP Dark Matter, and
Majorana Seesaw Mechanism

Xuewen Liu  Department of Physics, Yantai University, Yantai 264005, China
Shu-Yuan Guo  Department of Physics, Yantai University, Yantai 264005, China
Bin<EMAIL_ADDRESS>Department of Physics, Yantai University,
Yantai 264005, China Ying<EMAIL_ADDRESS>Department of Physics, Yantai
University, Yantai 264005, China

###### Abstract

We study a minimal extension of the Standard Model by introducing three right-
handed neutrinos and a new scotogenic scalar doublet, in which the mass
splittings between neutral and charged components are responsible for the
$W$-boson mass newly measured by the CDF collaboration. This model can not
only generate non-vanishing Majorana neutrino masses via the interaction of
right-handed neutrinos and scotogenic scalars, but also explain the Universe’s
missing matter in the form of FIMP dark matter. We also study the influence of
the mass splitting on the first order electroweak phase transition, and find
that it can further enhance the transition strength and thus induce
gravitational waves during the phase transition, which may be detected in the
forthcoming detectors such as U-DECIGO.

Keywords: $W$ boson mass, Dark matter, Neutrino, Electroweak phase transition,
Gravitational waves

##  1 Introduction

With high precision, the Standard Model (SM) explains the interactions of all
known fundamental particles. Despite of intensive investigations, no
significant deviations from the SM have been identified in the experiments,
suggesting that the SM constitutes the complete description of Nature.
However, several long-standing issues imply that new physics beyond the SM is
inevitable. The origin of neutrino masses and the formation of cold dark
matter are the two primary puzzles that any extensions of the SM should
address. Intriguingly, the newly measured $W$ boson mass provides yet another
impetus for new physics. Very recently, with the improved parton distribution
functions of the (anti)proton and new track reconstruction, the CDF
collaboration at Fermilab has released the world’s most precise direct
measurement of the $W$ boson mass [1], $m_{W}^{\mathrm{CDF}}=80.4435\pm
0.0094~{}\mathrm{GeV}$, based on $8.8~{}\mathrm{fb}^{-1}$ of data gathered
between $2002$ and $2011$, which is approximately $7\sigma$ away from the SM
prediction $m_{W}^{\rm EW}=80.3545\pm 0.0059$ GeV [2]. Furthermore, there is a
significant tension between the new CDF result and the direct measurements
implemented by the D0 collaboration at the Tevatron [3] and the ATLAS/LHCb
collaboration at the Large Hadron Collider (LHC) [4, 5] while the latter ones
are in good agreement with $m_{W}^{\rm EW}$. Although it is yet premature to
confirm the existence of new physics, a number of analyses on the new data and
interpretations in terms of physics beyond the SM have been carried out in the
literature [6, 7, 8, 9, 10, 11, 12, 13, 14].

In this work, we aim to bind the three issues together and interpret them in a
single setup, where we investigate the scotogenic model consisting of three
right-handed neutrinos and one $SU(2)$ doublet. This matter content was first
postulated in [15], which uses the scotogenic scalar as dark matter, and
explains the origin of tiny neutrino masses via a loop-induced process
involving the new matter contents. We found that the mass splitting between
the charged and neutral components of the $SU(2)$ doublet behaves like a
bridge to connect the three issues naturally. The scotogenic model is one of
the simplest loop realizations of the dim-5 Weinberg operator [16]. It lowers
the seesaw scale from GUT to the electroweak (EW) scale, meanwhile provides
candidate for dark matter in our Universe. The presence of a scotogenic scalar
also explains the latest $W$ boson mass discrepancy. This is mainly because
the authors in [7] demonstrated that the oblique parameters $S$, $T$, and $U$
must deviate from the SM estimate to establish a consistent EW global fit,
indicating the presence of mass-splitting in the scalar sector. Since the
scotogenic scalar has an $SU(2)$ representation, it yields a non-vanishing
contribution to oblique parameters.

The scalar sector in scotogenic model is also well known as the inert doublet
model [17, 18, 19, 20, 21] which naturally includes a DM candidate, i.e.
either the CP even Higgs or CP odd Higgs. However, the main difficulty in the
original inert model is to produce the correct relic density of dark matter
under the direct detection limits. Only in the Higgs funnel region (and heavy
mass region), can inert scalar annihilation attain the proper relic density.
On the other hand, the right-handed neutrino in this model can become a
promising dark matter candidate in most parameter space. The distinctive
structure of Yukawa couplings is not only helpful to achieve tiny neutrino
mass but beneficial to the freeze-in production of dark matter. The Feebly
Interacting Massive Particle (FIMP) couples to the thermal bath extremely
weak, so that it can not retain chemical equilibrium with the thermal bath.
Instead, the DM particles were produced by the decay processes from the scalar
sector, which provides the correct relic abundance. In addition, the new
introduced particles could also be used to explain the anomalies in flavor
physics [22].

The predictability of this paradigm is an essential part of its
attractiveness. Including the scotogenic scalar, the extra degrees of freedom
lead to an inevitable first order electroweak phase transition (EWPT) in the
early Universe. The scalar mass splittings required for explaining the $W$
boson mass happen to affect the magnitude of phase transition. EWPT in the
scotogenic scalar sector has been extensively studied [20, 21, 23, 24].
However, the new $W$ boson mass anomaly strongly motivates us to re-study the
brand new parameter space. As a consequence, the strong enough phase
transitions assure the generation of gravitational waves (GWs). Astonishingly,
we can use the resulting gravitational wave created to dig insight into the
$W$ bosons with the future facilities. The essential characteristics and
forecasts of the scotogenic model are presented in this work.

##  2 Neutrino mass and dark matter

The scotogenic model [15] is one of the simplest ways to link neutrino mass to
dark matter. Besides the particle content in SM, an additional doublet scalar
(denoted as $H_{2}$) and three generations of right-handed neutrinos
$N_{k}(k=1,2,3)$ are introduced. The new particles have odd parity under a new
$Z_{2}$ symmetry, while the SM particles are all $Z_{2}$ even. With $Z_{2}$
odd nature, the right-handed neutrinos cannot form Dirac masses
$\bar{L}\tilde{H}_{1}N+\rm{h.c.}$, where $H_{1}$ represents the SM Higgs
doublet and $\tilde{H}_{1}\equiv\epsilon H_{1}^{\ast}$. The relevant terms in
the Lagrangian, concerning right-handed neutrinos, are written as

$-\mathcal{L}_{N}=\frac{1}{2}m_{N}\overline{N^{c}}N+y_{N}\bar{L}\tilde{H}_{2}N+\rm{h.c.}.$
(1)

Similar to the type-I seesaw, the right-handed neutrinos could have Majorana
masses. The scalar potential, obeying SM gauge symmetries plus the discrete
$Z_{2}$, is given as

$\displaystyle
V=-\mu_{1}^{2}|H_{1}|^{2}+\mu_{2}^{2}|H_{2}|^{2}+\lambda_{1}|H_{1}|^{4}+\lambda_{2}|H_{2}|^{4}+\lambda_{3}|H_{1}|^{2}|H_{2}|^{2}+\lambda_{4}|H_{1}^{\dagger}H_{2}|^{2}+\frac{\lambda_{5}}{2}\left\\{(H_{1}^{\dagger}H_{2})^{2}+\rm{h.c.}\right\\}.$
(2)

The conserved $Z_{2}$ parity forbids $H_{2}$ from developing a vacuum
expectation value after EW symmetry breaking. However, we should mention that
the existence of an additional doublet alters the EWPT, which imprints the
detection signature as gravitational waves. The two doublets could be
expressed as,

$H_{1}=\left(\begin{array}[]{c}\phi_{1}^{+}\\\
\frac{1}{\sqrt{2}}(v+h+i\chi)\end{array}\right),H_{2}=\left(\begin{array}[]{c}H^{+}\\\
\frac{1}{\sqrt{2}}(H+iA)\end{array}\right),$ (3)

$h$ would be the only physical component of $H_{1}$, and it plays the role of
the observed Higgs scalar. For $H_{2}$, the four components are all physical,
with two neutral scalars, $H$(CP-even) and $A$(CP-odd), and two charged
scalars, $H^{\pm}$. Generally, the quartic couplings can be written in terms
of the physical scalars masses and $\mu_{2}$,

$\displaystyle\lambda_{3}$ $\displaystyle=$
$\displaystyle\frac{2}{v^{2}}\left(m_{H^{\pm}}^{2}-\mu_{2}^{2}\right),$ (4)
$\displaystyle\lambda_{4}$ $\displaystyle=$
$\displaystyle\frac{\left(m_{H}^{2}+m_{A}^{2}-2m_{H^{\pm}}^{2}\right)}{v^{2}},$
(5) $\displaystyle\lambda_{5}$ $\displaystyle=$
$\displaystyle\frac{\left(m_{H}^{2}-m_{A}^{2}\right)}{v^{2}}.$ (6)

We see that the $\lambda_{5}$ term controls the mass difference between the
two neutral components, and together with $\lambda_{4}$, the mass differences
in $H_{2}$ are determined. The couplings $\lambda_{i}$ are constrained by the
unitarity, vacuum stability, and perturbativity requirements (see [19] for
details), and in the following analysis we have taken these constraints into
account.

Now one could give a calculation on the neutrino mass in the scotogenic model,
which is given by [25]

$\displaystyle\left(\mathcal{M}_{\nu}\right)_{\alpha\beta}=\sum_{k=1}^{3}\frac{y_{N}^{\alpha
k}y_{N}^{\beta
k}}{32\pi^{2}}m_{N_{k}}\Big{[}\frac{m_{H}^{2}}{m_{H}^{2}-m_{N_{k}}^{2}}\ln\Big{(}\frac{m_{H}^{2}}{m_{N_{k}}^{2}}\Big{)}-\frac{m_{A}^{2}}{m_{A}^{2}-m_{N_{k}}^{2}}\ln\Big{(}\frac{m_{A}^{2}}{m_{N_{k}}^{2}}\Big{)}\Big{]},$
(7)

here $\alpha,\beta$ label the neutrino flavor indices. Thus the mass
difference between $H$ and $A$ plays an important role to produce the loop
neutrino mass, under the degenerate limit, i.e. $m_{H}=m_{A}$(or
$\lambda_{5}=0$), the loop mass would be vanished. If $\lambda_{5}\ll 1$, the
mass matrix could be rewritten as

$\displaystyle\left(\mathcal{M}_{\nu}\right)_{\alpha\beta}=\frac{\lambda_{5}v^{2}}{32\pi^{2}}\sum_{k=1}^{3}y_{N}^{\alpha
k}y_{N}^{\beta
k}\frac{m_{N_{k}}}{m_{0}^{2}-m_{N_{k}}^{2}}\times\Big{[}1-\frac{m_{N_{k}}^{2}}{m_{0}^{2}-m_{N_{k}}^{2}}\ln\Big{(}\frac{m_{0}^{2}}{m_{N_{k}}^{2}}\Big{)}\Big{]},$
(8)

with $m_{0}^{2}=(m_{H}^{2}+m_{A}^{2})/2$. The EW scale realization of neutrino
mass at the present energy frontier is relevant to the low-energy experiment.
Hence we focus on the masses of particles in loop below TeV. Apart from that,
we would pay our attention to thermal freeze-in production of the dark matter,
which is the lightest right-handed neutrino $N_{1}$. The out-of-equilibrium
accumulation, which we discuss in the dark matter section, requires relevant
coupling $y_{N}^{\alpha 1}$ to be extremely small. Hence $N_{1}$ contribution
to neutrino mass is negligible. Such a scenario is still consistent with the
neutrino oscillation result, since the converting probabilities care for the
squared mass differences of neutrinos, while one vanished neutrino mass is
allowed. One could estimate that to realize the neutrino mass at the TeV
scale, a combining constraint is $\lambda_{5}(y_{N}^{\alpha 2,3})^{2}\sim
10^{-11}$ [25]. It could be reached by tiny $\lambda_{5}$ while large
$y_{N}^{\alpha 2,3}$. But actually a too large Yukawa, e.g. $y_{N}>0.1$, would
result in strong contradiction with the non-observation of lepton flavor
violation(LFV) decays, e.g. $\ell_{i}\to\ell_{j}\gamma$.

The signature searching on collider of the new particles depends on relative
sizes of their masses. If the charged scalar is the next-to-lightest $Z_{2}$
odd particle, it has no choice but to decay into $N_{1}$, leaving a highly
ionizing charged track in the detector. The null result of searching for such
kind of track from heavy stable charged particle, on LHC, could exclude
$m_{H^{\pm}}$ up to $\sim 500~{}\rm{GeV}$ for $m_{N_{1}}\geq 1~{}\rm{MeV}$. If
the next-to-lightest $Z_{2}$ odd particle is $N_{2}$, the searching signature
would be two prompt charged leptons plus missing energy, which is from the
decay of $H^{\pm}\to N_{2}\ell^{\pm}$. The small decay rate of
$N_{2}\to\ell\bar{\ell}N_{1}$ would result in a decay-length which largely
exceeds size of the detector. Searching for two prompt leptons plus missing
energy on LHC could exclude $m_{H^{\pm}}$ up to $\sim 160~{}\rm{GeV}$, when
assuming the prompt charged leptons are either $e^{\pm}$ or $\mu^{\pm}$. The
exclusion limit would become weaker when decay into tau lepton is allowed. If
$N_{3}$ is lighter than the charged scalar, the decay of $H^{\pm}\to
N_{3}\ell^{\pm}$ would open and $N_{3}$ could subsequently decay into
$\ell\bar{\ell}N_{2}$. The relevant Yukawa couplings are strongly constrained
by the non-observation of $\mu\to e\gamma$, hence may result in a signature of
displaced lepton pair. The null result on LHC of the displaced lepton pair
could set constraints on masses of $N_{2,3}$ and $H^{\pm}$, depends on the
magnitude of Yukawa couplings. Decay products of the neutral scalars are
similar to the charged one except that charged leptons are replaced by
neutrinos, hence the limits on $m_{H/A}$ are weaker than the charged one. As
we focus on the case that charged scalar is the lightest $Z_{2}$ odd scalar,
adopting the above constraints for $H^{\pm}$ are enough for the scalar sector.
For detailed discussion of the collider searching, one is referred to [26]. In
the following discussions, we will adopt the corresponding constraints on
masses whenever the mass spectra are fixed. Besides, requirements on neutrino
mass and obeying LFV decays constraints are also imposed.

The discrete $Z_{2}$ parity on one hand gives loop masses for neutrinos, on
the other hand, it also suggests a candidate particle for dark matter. In
principle, both the lightest neutral scalar and sterile neutrino can be dark
matter. For the former situation, the neutral scalar could be thermally
freeze-out produced and becomes a typical WIMP. However, due to the current
stringent limits from the direct detection and relic density, most of the
parameter space has been ruled out of the WIMPs[24]. Another possibility is
that the lightest right-handed neutrino, i.e. $N_{1}$, acts as a FIMP dark
matter. In this situation, dark matter $N_{1}$ is produced out of thermal
equilibrium. Specifically speaking, $N_{1}$ is accumulated from the decay of
$X(H,A,H^{\pm})\to N_{1}L$, which relates to Yukawa $y_{N}^{\alpha 1}$
($y_{1}$ for short thereafter). The out of equilibrium condition requires the
decay rate to be smaller than the rate of Universe expansion, i.e.
$\Gamma(X\to N_{1}L)<H(T\sim m_{X})$, which would force $y_{1}$ to be small.
The dark matter $N_{1}$ yield, $Y_{N_{1}}(T)=n_{N_{1}}(T)/s(T)$, can be
computed by solving the following Boltzmann equation [27]

$sT\frac{dY_{N_{1}}}{dT}=-\frac{\gamma_{N_{1}}(T)}{H(T)},$ (9)

with $\gamma_{N_{1}}(T)$ representing the thermal averaged FIMP production
rate

$\gamma_{N_{1}}(T)=\sum_{X}\frac{g_{X}m_{X}^{2}T}{2\pi^{2}}K_{1}(m_{X}/T)\Gamma(X\to
N_{1}\ell),$ (10)

$g_{X}$ is the internal degrees of freedom of $X$, $s$ is the entropy density
of the Universe, $H(T)$ is the expansion rate of the Universe at a given
temperature and $K_{1}(x)$ is the Bessel function of the second kind.

The relevant study can be found in Ref. [25]. $N_{1}$ production will be
dominated by the decays of the scalars ($H,A,H^{\pm}$) while they are in
equilibrium with the thermal bath. The accumulation of $N_{1}$ could, in
principle, also from decay of heavier sterile neutrinos, i.e. $N_{2,3}\to
N_{1}\bar{\ell}\ell$, but Ref. [25] has verified these decays are subdominate.

The decay rates that enter into $\Gamma\left(X\to N_{1}\,L\right)$ are
calculated as,

$\displaystyle\Gamma(H\to N_{1}\bar{\nu})$
$\displaystyle=\frac{(m_{H}^{2}-m_{N_{1}}^{2})^{2}}{32\pi
m_{H}^{3}}y^{2}_{1},$ (11) $\displaystyle\Gamma(A\to N_{1}\,\bar{\nu})$
$\displaystyle=\frac{(m_{A}^{2}-m_{N_{1}}^{2})^{2}}{32\pi
m_{A}^{3}}y^{2}_{1},$ (12) $\displaystyle\Gamma(H^{+}\to N_{1}\,\bar{\ell})$
$\displaystyle=\frac{(m_{H^{+}}^{2}-m_{N_{1}}^{2})^{2}}{32\pi
m_{H^{+}}^{3}}y^{2}_{1}.$ (13)

By simply numerical calculations we can obtain the observed DM relic density
with the Yukawa coupling in the range $y_{1}\in[10^{-12},10^{-11}$], for tens
to hundreds GeV FIMPs. To study the correlations with the new CDF measured
$m_{W}$ and EWPT, we take $m_{H^{\pm}}=300$ GeV, $y_{1}=1.48\times 10^{-12}$
and $m_{N_{1}}=100$ GeV as a benchmark. Then we find the viable parameter
space for generating the relic density $\Omega h^{2}=0.12\pm 0.0012$ [28, 29],
which is shown by the purple band on mass splittings plane, see Fig. 1.

##  3 $W$ boson mass and Gravitational Waves

The relationship between $W$ boson mass and the oblique corrections is given
by [30],

$\displaystyle m_{W}^{2}=m_{W}^{2}(\mathrm{SM})+\frac{\alpha
c^{2}}{c^{2}-s^{2}}m_{Z}^{2}[-\frac{1}{2}\Delta S+c^{2}\Delta
T+\frac{c^{2}-s^{2}}{4s^{2}}\Delta U]\,,$ (14)

where $c=\cos\theta_{W}$, $s=\sin\theta_{W}$. To generate suitable $W$ boson
mass, we resort to a significant contribution from $\Delta S$ and $\Delta T$.
The expressions are [19]

$\displaystyle\Delta
S=\frac{1}{2\pi}\Big{[}\frac{1}{6}\log(\frac{m^{2}_{H}}{m^{2}_{H^{\pm}}})-\frac{5}{36}+\frac{m^{2}_{H}m^{2}_{A}}{3(m^{2}_{A}-m^{2}_{H})^{2}}+\frac{m^{4}_{A}(m^{2}_{A}-3m^{2}_{H})}{6(m^{2}_{A}-m^{2}_{H})^{3}}\log(\frac{m^{2}_{A}}{m^{2}_{H}})\Big{]},$
(15)

and

$\displaystyle\Delta T=\frac{1}{32\pi^{2}\alpha
v^{2}}\Big{[}F(m^{2}_{H^{\pm}},m^{2}_{A})+F(m^{2}_{H^{\pm}},m^{2}_{H})-F(m^{2}_{A},m^{2}_{H})\Big{]}$
(16)

where the function $F$ is defined by

$F(x,y)=\left\\{\begin{array}[]{lr}\frac{x+y}{2}-\frac{xy}{x-y}\log(\frac{x}{y}),&x\neq
y,\\\ 0,&\,x=y.\end{array}\right.$ (17)

The $T$ parameter vanishes when $m_{H^{\pm}}=m_{H}$ or $m_{H^{\pm}}=m_{A}$,
since these conditions lead to an exact custodial $SU(2)$ symmetry, in which
one of the neutral scalars joins the charged scalars to create a $SU(2)$
triplet. So to get proper oblique corrections, one knows that the key
characters are the mass splittings of the $Z_{2}$ odd scalars, namely $\Delta
M_{1}\equiv m_{H}-m_{H^{\pm}}$ and $\Delta M_{2}\equiv m_{A}-m_{H^{\pm}}$. In
previous studies of scalar multiplets models, the mass splittings are limited
by the previous $S,~{}T$ parameters fitted by the $m_{W}$ mass of PDG. Now we
have the right to reexamine the whole picture. In terms of $\Delta M_{1}$ and
$\Delta M_{2}$ we could directly get the CDF $W$ mass (the central value) by
using Eq. 14, which is shown in Fig. 1 by the black dash-dotted line in
$\Delta M_{1}$ and $\Delta M_{2}$ plane. We take $m_{H^{\pm}}=300$ GeV as a
benchmark model. In order to find the viable parameter regions, one should
also consider the correlations to the $W$ mass from the EW parameters. So we
use the global EW best fit values of $\Delta S$ and $\Delta T$ derived in Ref.
[7] (Table-III) to perform a fit by using $\chi^{2}({\bf O})=({\bf y}-\mu({\bf
O}))^{T}{\bf C}^{-1}({\bf y}-\mu({\bf O}))$, where ${\bf y}$ is the vector of
central values and ${\bf C}$ is the covariance matrix. In Fig. 1 we show 1-and
2-sigma regions allowed by $\Delta S$ and $\Delta T$ parameters in the mass
splittings plane. The $\Delta T$ parameter becomes negative if the charged
Higgs mass $m_{H^{\pm}}$ falls into the range between the masses of the two
neutral scalars, $m_{A}$ and $m_{H}$. As a result, either
$m_{H^{\pm}}>m_{A,H}$ or $m_{H^{\pm}}<m_{A,H}$ are permitted. Because of
$S\sim\log(m_{H,A}/m_{H^{\pm}})$, the $S$ parameter is disposed to be negative
in the former scenario, which is disfavored by the global fit [7]. Then we
only consider the case that charged Higgs is the lightest one in new higgs
doublet.

The rich spectrum of the scalar sector also drives the evolution dynamics of
the vacuum state non-trivial and provides the first order EWPT in the early
Universe. By considering the thermal loop effects we use the approximated
effective potential near the critical temperature [31]

$\displaystyle V_{\rm
eff}(h,T)\approx\frac{1}{2}(-\mu^{2}+cT^{2})h^{2}-\frac{\varepsilon
T}{12\pi}(h^{2})^{\frac{3}{2}}+\frac{\lambda}{4}h^{4},$ (18)

where the coefficient $\varepsilon$ quantifies the interactions between the
extra scalars and the Higgs boson. We have
$\varepsilon\approx(6m_{W}^{3}+3m_{Z}^{3})/{v^{3}}+2(\lambda_{3}/{2})^{3/2}+((\lambda_{3}+\lambda_{4}-\lambda_{5})/{2})^{3/2}+((\lambda_{3}+\lambda_{4}+\lambda_{5})/2)^{3/2}$
and
$c\approx(6m_{t}^{2}+6m_{W}^{2}+3m_{Z}^{2}+\frac{3}{2}m_{H}^{2})/(12v^{2})+(2\lambda_{3}+\lambda_{4})/12$.

The sizable oblique corrections favor large mass splittings $\Delta M_{1}$ and
$\Delta M_{2}$ in order to interpret the new measured $W$ boson mass, as shown
in Fig. 1. So the quartic couplings $\lambda_{4},\lambda_{5}$ in the effective
potential directly link the new $W$ deviation and the EWPT, then lead to an
affected EWPT.

The symmetry breaking phase starts to nucleate right after the Universe
cooling to the nucleation temperature $T_{n}$. The nucleation temperature is
defined by the equality of the nucleation rate per Hubble volume and the
Universe expansion rate, i.e. $\Gamma(T_{n})=H^{4}(T_{n})$, where
$\Gamma(T)\sim T^{4}e^{-S_{3}/T}$ is the decay rate per unit volume and
$S_{3}$ is the classical action for the $O(3)$ symmetric bounce solution [32].
For a radiation-dominated Universe and a first order phase transition (FOPT)
happening around the EW scale, $T_{n}$ can be solved by $S_{3}/T_{n}\sim 140$
[33]. This criterion will be taken as the sufficient condition for a FOPT in
this work. Further we use the package cosmoTransition [34] to calculate the
bounce solution and $T_{n}$ for the effective potential $V_{\rm eff}$. The
custodial symmetry in the scalar potential leads to a mass degeneration of
$m_{A}$ and $m_{H^{\pm}}$ in previous research (or before the $W$ mass
deviation is identified). Generally, the EWPT strength depends on the mass
splitting since it affects the potential barrier from custodial symmetry
breaking. Now we have the chance to open the new parameter space by taking
$m_{A}\neq m_{H^{\pm}}$, we can find that such mass splitting $\Delta M_{2}$
enhances the EWPT strength in an obvious way, as shown in Fig. 1 with red-
dashed contours. The reason is mainly that the larger $\Delta M_{1}$ and/or
$\Delta M_{2}$ leads to a higher barrier between the symmetric vacuum and the
symmetry breaking one. (see Eqs. 5, 6 and Eq. 18). This conclusion could also
be generalized to other scalar multiplets models.

Figure 1: (Color online) Best fit region for $\Delta S,\Delta T$ operators,
FOPT strength $v_{n}/T_{n}$ and DM relic density, as functions of $\Delta
M_{1}$ and $\Delta M_{2}$ where $\Delta M_{1}\equiv m_{H}-m_{H^{\pm}}$ and
$\Delta M_{2}\equiv m_{A}-m_{H^{\pm}}$. The green and yellow bands correspond
to the best fit with 1- and 2-sigma regions required by the new CDF measured
$m_{W}$. The black dash-dotted line satisfies the central value of the CDF $W$
mass. In the red region the strongly FOPT can be realized and the red dashed
lines represent the contours of $v_{n}/T_{n}$. The gray line represents the
condition $m_{H}=m_{A}$ which means that one cannot obtain the neutrino mass
here. While the neutrino mass can be generated in the most parameter space
with adjusted Yukawa couplings. In the purple region, one can yield the
correct DM relic density $\Omega h^{2}=0.12\pm 0.0012$ [28], with coupling
$y_{1}=1.48\times 10^{-12}$ and $m_{N_{1}}=100$ GeV. The red point is the
selected benchmark point (BP) to investigate GWs physics below, which is:
$m_{H}=423$ GeV, $m_{A}=426$ GeV, $m_{H^{\pm}}=300$ GeV, $\lambda_{3}=1.5$,
$\lambda_{4}=2.98$, $\lambda_{5}=-0.04$, $T_{n}=158.78$ GeV,
$v_{n}/T_{n}=4.1$.

During a FOPT, stochastic GWs come from three sources: bubble collisions [35,
36], sound waves in the plasma [37, 38] and the magneto-hydrodynamics
turbulence [39, 40]. Recent studies show that the bubble collision
contribution to the GWs can generally be neglected because only a tiny
fraction of the FOPT energy deposits in the bubble wall[41]. It turns out that
the dominant contribution comes from the sound waves as most FOPT energy is
pumped into the surrounding fluid shells [42]. Then the turbulence is another
main source after a finite period of the sound wave[43]. Consequently, the GWs
spectrum today can be expressed as (see details in supplementary file.)

$\Omega_{\rm GW}(f)=\Omega_{\rm sw}(f)+\Omega_{\rm turb}(f),$ (19)

where $f$ is the frequency, the subscripts “sw” and “turb” denote sound waves
and turbulence respectively.

The GWs produced from a FOPT around EW scale have the potential to be probed
by the next generation space-based laser interferometers such as LISA [44],
BBO [45], TianQin [46, 47], Taiji [48, 49, 50], DECIGO [51, 52] and U-DECIGO
[53]. In Fig. 2 we plot the GW spectrum with the benchmark point. We find that
the GW signals are within the detectability of U-DECIGO, and this possibility
gives a complementary detection to the other searches. While to produce
stronger signals one may consider the extension of the scalar sector, which
further enhances the sensitivities [54].

Figure 2: (Color online) GW spectrum $\Omega_{\rm GW}h^{2}$ from the benchmark
point where $\alpha=0.021$, ${\beta}/{H}=722.9$, depicted with the red line.
The colored regions represent the expected sensitivities of the projected GWs
interferometers.

##  4 Conclusion

Eventually, we successfully explain our Universe with the new measured $W$
boson mass, the correct abundance of FIMP dark matter, and the Majorana
neutrino masses within a minimal extension of the SM. We found that both
neutrino mass and $W$ boson mass are sensitive to the non-vanishing $\Delta
M_{1}$ and $\Delta M_{2}$, which also influences the freeze-in production of
sterile neutrino dark matter. As a result, the non-degeneration between the
charged and neutral components of a scotogenic scalar naturally links the
three fundamental issues. More importantly, this model makes a precise
prediction that will be validated in future GW experiments. The additional
multiplet in the scalar sector drives the evolution dynamics of the vacuum
state non-trivial. Thus in the early Universe, it provides a first-order
electroweak phase transition, which is further enhanced by the new required
mass splitting $m_{A}-m_{H^{\pm}}$. This is thought to be an efficient
approach to generating detectable GWs. The future gravitational wave detectors
will help reveal valuable information about the nature of this model.

## Acknowledgements

This work was supported by the National Natural Science Foundation of China
under the grants Nos. 11805161, 12005180, and 11975195, by the Natural Science
Foundation of Shandong Province under Grant No. ZR2020QA083 and ZR2019JQ04 and
by the Project of Shandong Province Higher Educational Science and Technology
Program under Grants No. 2019KJJ007.

## References

  * [1] CDF Collaboration, T. Aaltonen et al., High-precision measurement of the $W$ boson mass with the CDF II detector, Science 376 (2022), no. 6589 170–176.
  * [2] J. de Blas, M. Ciuchini, E. Franco, A. Goncalves, S. Mishima, M. Pierini, L. Reina, and L. Silvestrini, Global analysis of electroweak data in the Standard Model, [arXiv:2112.07274].
  * [3] D0 Collaboration, V. M. Abazov et al., Measurement of the W Boson Mass with the D0 Detector, Phys. Rev. Lett. 108 (2012) 151804, [arXiv:1203.0293].
  * [4] ATLAS Collaboration, M. Aaboud et al., Measurement of the $W$-boson mass in pp collisions at $\sqrt{s}=7$ TeV with the ATLAS detector, Eur. Phys. J. C 78 (2018), no. 2 110, [arXiv:1701.07240]. [Erratum: Eur.Phys.J.C 78, 898 (2018)].
  * [5] LHCb Collaboration, R. Aaij et al., Measurement of the W boson mass, JHEP 01 (2022) 036, [arXiv:2109.01113].
  * [6] Y.-Z. Fan, T.-P. Tang, Y.-L. S. Tsai, and L. Wu, Inert Higgs Dark Matter for New CDF W-boson Mass and Detection Prospects, [arXiv:2204.03693].
  * [7] C.-T. Lu, L. Wu, Y. Wu, and B. Zhu, Electroweak Precision Fit and New Physics in light of $W$ Boson Mass, [arXiv:2204.03796].
  * [8] P. Athron, A. Fowlie, C.-T. Lu, L. Wu, Y. Wu, and B. Zhu, The $W$ boson Mass and Muon $g-2$: Hadronic Uncertainties or New Physics?, [arXiv:2204.03996].
  * [9] G.-W. Yuan, L. Zu, L. Feng, and Y.-F. Cai, $W$-boson mass anomaly: probing the models of axion-like particle, dark photon and Chameleon dark energy, [arXiv:2204.04183].
  * [10] A. Strumia, Interpreting electroweak precision data including the $W$-mass CDF anomaly, [arXiv:2204.04191].
  * [11] J. M. Yang and Y. Zhang, Low energy SUSY confronted with new measurements of W-boson mass and muon g-2, [arXiv:2204.04202].
  * [12] J. de Blas, M. Pierini, L. Reina, and L. Silvestrini, Impact of the recent measurements of the top-quark and W-boson masses on electroweak precision fits, [arXiv:2204.04204].
  * [13] C.-R. Zhu, M.-Y. Cui, Z.-Q. Xia, Z.-H. Yu, X. Huang, Q. Yuan, and Y. Z. Fan, GeV antiproton/gamma-ray excesses and the $W$-boson mass anomaly: three faces of $\sim 60-70$ GeV dark matter particle?, [arXiv:2204.03767].
  * [14] A. Addazi, A. Marciano, A. P. Morais, R. Pasechnik, and H. Yang, CDF II $W$-mass anomaly faces first-order electroweak phase transition, [arXiv:2204.10315].
  * [15] E. Ma, Verifiable radiative seesaw mechanism of neutrino mass and dark matter, Phys. Rev. D 73 (2006) 077301, [hep-ph/0601225].
  * [16] S. Weinberg, Baryon and Lepton Nonconserving Processes, Phys. Rev. Lett. 43 (1979) 1566–1570.
  * [17] N. G. Deshpande and E. Ma, Pattern of Symmetry Breaking with Two Higgs Doublets, Phys. Rev. D 18 (1978) 2574.
  * [18] L. Lopez Honorez and C. E. Yaguna, A new viable region of the inert doublet model, JCAP 01 (2011) 002, [arXiv:1011.1411].
  * [19] A. Arhrib, R. Benbrik, and N. Gaur, $H\to\gamma\gamma$ in Inert Higgs Doublet Model, Phys. Rev. D 85 (2012) 095021, [arXiv:1201.2644].
  * [20] D. Borah and J. M. Cline, Inert Doublet Dark Matter with Strong Electroweak Phase Transition, Phys. Rev. D 86 (2012) 055001, [arXiv:1204.4722].
  * [21] G. Gil, P. Chankowski, and M. Krawczyk, Inert Dark Matter and Strong Electroweak Phase Transition, Phys. Lett. B 717 (2012) 396–402, [arXiv:1207.0084].
  * [22] Y. Li and C.-D. Lü, Recent Anomalies in B Physics, Sci. Bull. 63 (2018) 267–269, [arXiv:1808.02990].
  * [23] N. Blinov, S. Profumo, and T. Stefaniak, The Electroweak Phase Transition in the Inert Doublet Model, JCAP 07 (2015) 028, [arXiv:1504.05949].
  * [24] S. Fabian, F. Goertz, and Y. Jiang, Dark matter and nature of electroweak phase transition with an inert doublet, JCAP 09 (2021) 011, [arXiv:2012.12847].
  * [25] E. Molinaro, C. E. Yaguna, and O. Zapata, FIMP realization of the scotogenic model, JCAP 07 (2014) 015, [arXiv:1405.1259].
  * [26] A. G. Hessler, A. Ibarra, E. Molinaro, and S. Vogl, Probing the scotogenic FIMP at the LHC, JHEP 01 (2017) 100, [arXiv:1611.09540].
  * [27] L. J. Hall, K. Jedamzik, J. March-Russell, and S. M. West, Freeze-In Production of FIMP Dark Matter, JHEP 03 (2010) 080, [arXiv:0911.1120].
  * [28] Planck Collaboration, N. Aghanim et al., Planck 2018 results. VI. Cosmological parameters, Astron. Astrophys. 641 (2020) A6, [arXiv:1807.06209]. [Erratum: Astron.Astrophys. 652, C4 (2021)].
  * [29] M. Abdughani, Y.-Z. Fan, L. Feng, Y.-L. S. Tsai, L. Wu, and Q. Yuan, A common origin of muon g-2 anomaly, Galaxy Center GeV excess and AMS-02 anti-proton excess in the NMSSM, Sci. Bull. 66 (2021) 2170–2174, [arXiv:2104.03274].
  * [30] M. E. Peskin and T. Takeuchi, Estimation of oblique electroweak corrections, Phys. Rev. D 46 (1992) 381–409.
  * [31] D. J. H. Chung, A. J. Long, and L.-T. Wang, 125 GeV Higgs boson and electroweak phase transition model classes, Phys. Rev. D 87 (2013), no. 2 023509, [arXiv:1209.1819].
  * [32] A. D. Linde, Decay of the False Vacuum at Finite Temperature, Nucl. Phys. B 216 (1983) 421. [Erratum: Nucl.Phys.B 223, 544 (1983)].
  * [33] M. Quiros, Finite temperature field theory and phase transitions, pp. 187–259, 1, 1999.  hep-ph/9901312.
  * [34] C. L. Wainwright, CosmoTransitions: Computing Cosmological Phase Transition Temperatures and Bubble Profiles with Multiple Fields, Comput. Phys. Commun. 183 (2012) 2006–2013, [arXiv:1109.4189].
  * [35] S. J. Huber and T. Konstandin, Gravitational Wave Production by Collisions: More Bubbles, JCAP 09 (2008) 022, [arXiv:0806.1828].
  * [36] Y. Di, J. Wang, R. Zhou, L. Bian, R.-G. Cai, and J. Liu, Magnetic Field and Gravitational Waves from the First-Order Phase Transition, Phys. Rev. Lett. 126 (2021), no. 25 251102, [arXiv:2012.15625].
  * [37] 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 (2015), no. 12 123009, [arXiv:1504.03291].
  * [38] 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 (2021) 001, [arXiv:2007.08537].
  * [39] P. Binetruy, A. Bohe, C. Caprini, and J.-F. Dufaux, Cosmological Backgrounds of Gravitational Waves and eLISA/NGO: Phase Transitions, Cosmic Strings and Other Sources, JCAP 06 (2012) 027, [arXiv:1201.0983].
  * [40] C. Caprini, R. Durrer, and G. Servant, The stochastic gravitational wave background from turbulence and magnetic fields generated by a first-order phase transition, JCAP 12 (2009) 024, [arXiv:0909.0622].
  * [41] D. Bodeker and G. D. Moore, Electroweak Bubble Wall Speed Limit, JCAP 05 (2017) 025, [arXiv:1703.08215].
  * [42] J. Ellis, M. Lewicki, and J. M. No, On the Maximal Strength of a First-Order Electroweak Phase Transition and its Gravitational Wave Signal, JCAP 04 (2019) 003, [arXiv:1809.08242].
  * [43] J. Ellis, M. Lewicki, and J. M. No, Gravitational waves from first-order cosmological phase transitions: lifetime of the sound wave source, JCAP 07 (2020) 050, [arXiv:2003.07360].
  * [44] LISA Collaboration, P. Amaro-Seoane et al., Laser Interferometer Space Antenna, [arXiv:1702.00786].
  * [45] J. Crowder and N. J. Cornish, Beyond LISA: Exploring future gravitational wave missions, Phys. Rev. D 72 (2005) 083005, [gr-qc/0506015].
  * [46] TianQin Collaboration, J. Luo et al., TianQin: a space-borne gravitational wave detector, Class. Quant. Grav. 33 (2016), no. 3 035010, [arXiv:1512.02076].
  * [47] Y.-M. Hu, J. Mei, and J. Luo, Science prospects for space-borne gravitational-wave missions, Natl. Sci. Rev. 4 (2017), no. 5 683–684.
  * [48] W.-R. Hu and Y.-L. Wu, The Taiji Program in Space for gravitational wave physics and the nature of gravity, Natl. Sci. Rev. 4 (2017), no. 5 685–686.
  * [49] W.-H. Ruan, Z.-K. Guo, R.-G. Cai, and Y.-Z. Zhang, Taiji program: Gravitational-wave sources, Int. J. Mod. Phys. A 35 (2020), no. 17 2050075, [arXiv:1807.09495].
  * [50] Z.-W. Zhao, L.-F. Wang, J.-F. Zhang, and X. Zhang, Prospects for improving cosmological parameter estimation with gravitational-wave standard sirens from Taiji, Sci. Bull. 65 (2020), no. 16 1340–1348, [arXiv:1912.11629].
  * [51] S. Kawamura et al., The Japanese space gravitational wave antenna: DECIGO, Class. Quant. Grav. 28 (2011) 094011.
  * [52] S. Kawamura et al., The Japanese space gravitational wave antenna DECIGO, Class. Quant. Grav. 23 (2006) S125–S132.
  * [53] H. Kudoh, A. Taruya, T. Hiramatsu, and Y. Himemoto, Detecting a gravitational-wave background with next-generation space interferometers, Phys. Rev. D 73 (2006) 064006, [gr-qc/0511145].
  * [54] L. Bian and X. Liu, Two-step strongly first-order electroweak phase transition modified FIMP dark matter, gravitational wave signals, and the neutrino mass, Phys. Rev. D 99 (2019), no. 5 055003, [arXiv:1811.03279].