# Study of radio transients from the quiet Sun during an extremely quiet time

Surajit Mondal Center for Solar Terrestrial Research, New Jersey Institute of
Technology, Newark, New Jersey, 07102, United States of America  Divya Oberoi
National Centre for Radio Astrophysics, Tata Institute of Fundamental
Research, Pune-411007, India  Ayan Biswas National Centre for Radio
Astrophysics, Tata Institute of Fundamental Research, Pune-411007, India

(Received TBD ; Revised ; Accepted TBD )

###### Abstract

In this work we study a class of recently discovered metrewave solar
transients referred to as Weak Impulsive Narrowband Quiet Sun Emission
(WINQSEs, Mondal et al., 2020). Their strength is a few percent of the quiet
Sun background and are characterised by their very impulsive, narrow-band and
ubiquitous presence in quiet Sun regions. Mondal et al. (2020) hypothesised
that these emissions might be the radio counterparts of the nanoflares and
their potential significance warrants detailed studies. Here we present an
analysis of data from an extremely quiet time and with an improved methodology
over the previous work. As before, we detect numerous WINQSEs, which we have
used for their further characterisation. Their key properties, namely, their
impulsive nature and ubiquitous presence on the quiet Sun are observed in
these data as well. Interestingly, we also find some of the observed
properties to differ significantly from the earlier work. With this
demonstration of routine detection of WINQSEs, we hope to engender interest in
the larger community to build a deeper understanding of WINQSEs.

Quiet solar corona (1992); Solar corona (1483); Solar coronal radio emission
(1993); Solar coronal heating (1989)

††facilities: Murchison Widefield Array (MWA) (Lonsdale et al., 2009; Tingay
et al., 2013)††software: astropy (Astropy Collaboration et al., 2013, 2018),
matplotlib (Hunter, 2007), Numpy (Harris et al., 2020), SciPy (Virtanen et
al., 2020), CASA (McMullin et al., 2007)

##  1 Introduction

Radio transients from the Sun were discovered at very the dawn of solar radio
astronomy (e.g. Wild et al., 1963; Wild & Smerd, 1972, etc.). Since then
numerous studies have analysed and studied their various properties ranging
from emission mechanisms, spectro-temporal characteristics, relationships with
other solar phenomenon, signatures in other wavebands and so on (e.g.
Takakura, 1967; Sturrock, 1964; Melrose, 1980; McLean & Labrum, 1985; Pick &
Vilmer, 2008; Kontar et al., 2017; Reid & Kontar, 2017; Reid, 2020, and many
others). Most of these transients, particularly those happening in the meter
wavelengths regime, are associated with active regions and violent solar
activity like X-ray flares and coronal mass ejections (e.g. Pick & Vilmer,
2008; Bain & Fletcher, 2009; Krucker et al., 2011; Reid, 2020; Mondal &
Oberoi, 2021, etc.). However, there are hypotheses which suggest that such
radio transients, albeit much weaker, should also be present in the quiet
solar corona (Che, 2018). The primary reasoning behind this expectation comes
from the Parker’s nanoflare hypothesis (Parker, 1988) which states that the
solar corona maintains its million K temperature because of the energy
continuously deposited by numerous tiny nanoflares happening everywhere in the
corona. If this hypothesis is true, it implies that small scale reconnections
are happening everywhere and all the time in the corona. These small scale
reconnections, like their larger counterparts, will also accelerate electrons,
which in turn, on their interaction with the thermal plasma, should emit
plasma emission, observable in the radio bands (Che, 2018). To the best of our
knowledge, the first detection of the ubiquitous presence of these radio
transients from the quiet sun came from Mondal et al. (2020) (henceforth
referred to as M20). They found that these transients are highly impulsive
(duration $\lesssim 1$s) and narrowband in nature. Based on these properties
they also concluded that these emissions are generated by nonthermal
processes. They also found these emissions to be present throughout the quiet
sun. Mondal (2021) has demonstrated, although in a single instance, that the
energy associated with a group of these transients is $\sim 10^{25}$ergs.
These observations suggest the possibility that these transients are the radio
counterparts of the long hypothesised nanoflares. Henceforth we refer to these
emissions as Weak Impulsive Narrowband Quiet Sun Emissions (WINQSEs).

Due to the large potential significance of this discovery, it is crucial to
verify the presence of WINQSEs in other datasets obtained under different
solar conditions, using independent detection techniques and especially during
quiet Sun times. As a part of this exercise, recently Sharma et al. (2022)
(henceforth referred to as S22) have reported detection of WINQSEs using a
completely independent technique and analysis pipeline using a different
dataset. Interestingly S22 also identified a very weak type III solar radio
burst like drifting emission feature, which supports the hypothesis that
WINQSEs are the weaker cousins of the stronger type III bursts. Continuing
with the same motivation, here we present an analysis of data from an
extremely quiet time on 20th June, 2020. The solar activity levels on this day
are much lower than those prevailing during the observations of M20 and S22.
We also take this opportunity to improve the analysis methodology beyond what
was used in M20.

This paper is structured as follows. In Section 2 we present the observation
details, Section 3 discusses the details of data analysis including the
improvements over M20. In Section 4 we present the results, in Section 5 the
discussion and finally Section 6 gives the concluding remarks.

##  2 Observation details

The data presented here come from a solar observation with the Murchison
Widefield Array Phase II (MWA, Lonsdale et al., 2009; Tingay et al., 2013;
Wayth et al., 2018). These data were acquired on 20 June, 2020 between
03:35:00–04:40:00 UT under the project ID G0002. The MWA was in its extended
configuration, with a maximum baseline of $\sim 5$km. The observing frequency
spanned the band from 119.68–150.40 MHz. The data was acquired with a time and
frequency resolution of 0.5 s and 10 kHz respectively. Here we present the
analysis at four selected spectral slices, each of width 160 kHz. These slices
are centered around 120.22, 127.9, 135.58 and 139.9 MHz.

The Sun was extremely quiet on this
day111https://solarmonitor.org/index.php?date=20200620. No X-ray flares were
reported by the Geostationary Operational Environmental Satellites (GOES) and
no events of any kind were reported by the Space Weather Prediction Center
event list for the day. The Learmonth Spectrograph operating between 25–180
MHz also did not report any radio flare on this day. No NOAA active region was
present on the visible solar disc.

##  3 Data analysis

Making the large number of images needed for this investigation necessarily
requires a robust automated pipeline. Such a tool has recently been developed
and christened the Automated Imaging Routine for Compact Arrays for Radio Sun
(AIRCARS; Mondal et al., 2019). The credentials of this self-calibration based
pipeline have already been established in multiple recent works requiring high
quality spectroscopic snapshot imaging (e.g. Mohan et al., 2019; Mondal et
al., 2020; Mohan, 2021; Mondal & Oberoi, 2021, etc.). It is known from past
experience that the presence of strong compact sources on the Sun makes
calibration easier and the data corresponding to a large but featureless quiet
Sun is harder to calibrate (Mondal et al., 2019). While the imaging of the
data used here was done using AIRCARS, to improve the accuracy of calibration
for these quiet Sun data, an improved strategy from what was followed in M20
was employed. For this work, the calibration solutions were determined using
data averaged over 9 s. This time integration helps reduce the thermal noise
substantially and is especially useful in improving the signal-to-noise of
visibilities from longer baselines. The calibration solutions obtained were
linearly interpolated in time and applied to the 0.5 s resolution data which
were then imaged. The analysis for each of the four frequencies presented here
was carried out independently following this procedure. This approach relies
on the assumption that the solar variability over the integration period is
small enough to be ignored. We find this to be a good assumption for this
extremely quiet time and is evidenced by the higher imaging dynamic ranges of
the images produced using solutions derived from time-averaged visibilities.

On its own, AIRCARS does not ensure a common flux density scale across
independent observations. For the present analysis, this can lead to
discontinuities in the fluxscale at the boundaries of observations, typically
of a duration of 4-5 minutes, though this issue has now been resolved in its
next incarnation, named, Polarimetry using Automated Imaging Routine for
Compact Arrays for the Radio Sun (P-AIRCARS; Kansabanik et al., 2022b). A
precise absolute flux density is not required for this work and we have
adopted the strategy described below for relative flux density calibration
between all of the observations at a given frequency. It also takes care of
the amplitude variability of the order of 4-5% of the MWA instrumental
response which otherwise is not possible to take care of by using the method
described in Kansabanik et al. (2022a). While under most circumstances, this
level of variability is too small to be of any consequence, in the present
context it is important to correct for it as it is of a magnitude comparable
to the strengths of WINQSEs.

A correction for the instrumental primary beam is first applied to the solar
maps and then it is enforced that the median disc integrated flux density of
the Sun be the same for all five-minute chunks of data used. The underlying
assumption is that as the Sun was very quiet during our observing period, its
disc-integrated flux density should also be stable in this period. The
characteristics of WINQSEs - flux densities of $\sim$ mSFU (as compared to a
few SFU from the Sun) and short durations ($\lesssim$1 s) - imply that they
are too weak to influence the disc-integrated flux density in any significant
manner. We also found that the point-spread-function (psf) corresponding to
different time slices for a given frequency can change slightly across time.
Flux densities in radio maps are measured in units of $Jy/beam$, where $beam$
refers to the area of the psf. Hence, slight changes in the psf area can give
rise to low-level jumps in flux density estimates. To avoid contamination due
to this, we have smoothed all of the images used here to a common coarser
resolution of $280^{"}\times 280^{"}$.

A key step for identifying WINQSEs in M20 was to calculate the median image,
which requires accurate alignment of images relative to each other. M20
achieved this by aligning the radio peak of the active region present in their
dataset across images. Absence of a clearly identifiable compact feature
forces us to look for other alternatives to achieve this objective. We compute
a center of intensity ($\hat{X}$, $\hat{Y}$) for each 0.5 s solar image,
defined as:

$\hat{X}=\frac{\sum_{i}x_{i}}{\sum_{i}1},\quad\hat{Y}=\frac{\sum_{i}y_{i}}{\sum_{i}1},$
(1)

where ($x_{i}$, $y_{i}$) are the coordinates of the pixels which exceed the
threshold chosen visually such that no noise contours appear in the image and
the region it demarcates is reasonably symmetric. The centers of intensity for
all the images were aligned to lie at the same pixel. The aligned maps were
used to generate the median map, one for each five-minute data chunk and
frequency. An example 0.5 s map at 127.9 MHz and the corresponding median map
are shown in Fig. 1.

Figure 1: Left panel: An example 0.5 s map at 127.9 MHz. Right panel: The
corresponding median map at 127.9 MHz.

##  4 Results

For detection of WINQSEs in the large number of images used in this work, we
have modified the technique used in M20. First, all 0.5 s solar images were
smoothed using a circular kernel of size 280". Let the images prior and post
this smoothing operation be denoted by $I_{init}$ and $I_{flux}$ respectively.
As mentioned in the previous section, this smoothing operation ensured that
the numerical value of each pixel in the radio image $I_{flux}$ is equal to
the flux density for the resolution element corresponding to the circular
smoothing kernel centered on that pixel in $I_{init}$. The median map obtained
using $I_{init}$ images was used to obtain $I_{median,flux}$ using the same
technique as that used to obtain $I_{flux}$. $I_{median,flux}$ was then
subtracted from each $I_{flux}$ to obtain the residual maps for each 0.5 s
timeslice, henceforth referred to as $I_{res}$. The pixel values in $I_{res}$
is equivalent to $\Delta F_{i,\nu,t}=F_{i,\nu,t}-<F_{i,\nu}>$, where
$F_{i,\nu,t}$ is the flux density of a psf sized region $i$ at frequency $\nu$
some time $t$ and $<F_{i,\nu}>$ is the median flux density for region $i$ at
frequency $\nu$. A WINQSE is said to be present when the pixel value in
$I_{res}$ exceeds a given threshold. In this work, the threshold is set to
$3\sigma$, where $\sigma$ is the rms in the corresponding $I_{flux}$ in a
region away from the Sun. To ensure that only independent pixels are searched
for WINQSEs, the search algorithm steps through the pixels in steps of 7
pixels in both directions which translates to $280^{"}$ in each direction and
is equal to the effective image resolution. Figure 2 shows the contours of an
example residual image with a bright WINQSE, overlaid on a AIA 193Å image. The
compact feature seen in the image is the WINQSE detected in this snapshot and
it is evident that it has been detected with good signal-to-noise. The
following subsections discuss the various properties of the WINQSEs detected
in these data.

Figure 2: Radio contours of an example WINQSE at 136 MHz are overlaid on a 193
Å image from AIA/SDO. The contour levels are at 0.7, 0.85 and 0.99 times the
peak of the median subtracted map.

###  4.1 Flux density distribution

For convenience and consistency with M20, for a given frequency $\nu$, we
refer to the collection of $\Delta F_{i,\nu,t}$ for all psf sized regions $i$
s and all times $t$ s as $\Delta F$ and similarly to $<F_{i,\nu}>$ as $F$. We
define $\Delta F/F$ to be the ratio of the flux density in residual images
($I_{res}$) to that in the corresponding median image ($I_{median,flux}$) at
the same pixel location. Figure 3 shows the histogram of $\Delta F/F$ in
detected WINQSEs. These distributions are strikingly different from those
found by M20. In Sec. 5.2 we try to understand and reconcile these
differences.

The red line in Fig. 3 shows a lognormal fit to the data shown with blue
points. It is evident that the lognormal distribution is a good fit to the
data. However, it is worth noting that the part of the distribution where
$\Delta F/F\lesssim 0.1$ is affected by incompleteness due to a combination of
multiple factors including sensitivity, and spatial and temporal resolution.
These factors are not expected to significantly affect the data at higher
$\Delta F/F$ values and hence this part of the histogram is robust. The fact
that such lognormal distributions have been observed in earlier works when
studying weak impulsive emissions, though at EUV wavelengths and with rather
different temporal and spatial sampling (Pauluhn & Solanki, 2007; Tajfirouze &
Safari, 2012; Upendran & Tripathi, 2021), suggests that WINQSEs too might
follow a lognormal distribution. Testing this hypothesis in a more robust
manner requires an instrument with better sensitivity and resolution. We plan
to explore this with the next phase of the MWA (Phase III) which will offer a
factor of two larger effective collecting area.

Figure 3: Histograms of $\Delta F/F$ are shown using the blue circles. The red
line shows the lognormal fit to these data. The observation frequencies are
mentioned in each panel.

###  4.2 Temporal width distribution

The temporal width of a WINQSE is defined as the time duration for which the
value of a pixel in $I_{res}$ continuously exceeds the 3$\sigma$ threshold.
Figure 4 shows the temporal width distribution of the detected WINQSEs. The
powerlaw index of the distributions is close to -2, very similar to that found
in M20. The bursts are highly impulsive with a steep powerlaw like
distribution. The peak of the distribution lies at the instrumental resolution
of 0.5 s, implying that most of the WINQSEs are unresolved in time.

Figure 4: Temporal width distribution of the detected WINQSEs.

###  4.3 Wait time distribution

The wait time between WINQSEs is defined as the time interval between the
successive times when the pixel value at a given pixel in $I_{res}$ exceeds
the 3$\sigma$ threshold. The wait time distribution of detected WINQSEs is
shown in Fig. 5. The wait time distribution is very similar to that found in
M20. During the interpretation of the wait time distribution, a key point was
however missed by M20. Kivelä & Porter (2015) pointed out that the sharp
cutoff at high wait times, like what is observed in Fig. 5, can arise because
of data gaps and hence is not reliable. To take this into account, in Fig. 5,
we draw red and black dashed lines to indicate the 10% and 20% error
probabilities respectively. This implies that the data is sufficient to say
that at low wait times the distribution follows a powerlaw, but the nature of
the distribution at higher wait times cannot be constrained with these data.
The powerlaw nature of wait time distribution is observed for solar flares and
has been explained using a non-stationary Poisson model for solar flares
(Aschwanden & McTiernan, 2010). It is possible that similar behaviour is also
present for WINQSEs. M20 mention this aspect and alluded to the possibility
that the observed wait time distribution might be severely affected by a lack
of data. Based on the analysis here we conclude that this is indeed the case
for the high wait time regime.

Figure 5: Wait time distribution of the detected WINQSEs. The red and black
lines show the wait time till which the uncertainties due to incomplete
sampling is less than 10% and 20% respectively.

###  4.4 Spatial distribution

A key result of M20 was that the WINQSEs are ubiquitous in nature. Figure 6
shows the fractional occupancy of the detected WINQSEs over the radio Sun. The
fractional occupancy is defined as the fraction of observation time for which
WINQSEs were detected in a resolution element. There is no patch on the Sun
where WINQSEs are absent.

Figure 6: Fractional occupancy distribution of the detected WINQSEs. Contours
of the median map at respective frequencies has been overlaid. The contour
levels correspond to 0.2, 0.4, 0.6, 0.8, 0.9, 0.95 times the peak in the
median map.

##  5 Discussion

###  5.1 Improving the data analysis technique over M20

We have taken this opportunity to improve the analysis methodology to account
for two low-level effects which were realized during the course of this work
and are listed below:

  1. 1.

In M20 each five minute datachunk was processed independently, this can lead
to small differences in absolute fluxscale between data chunks. At a few
percent level, these differences are smaller than the uncertainty in the
absolute flux density and usually too small to be of any consequence. In the
present context, however, it is important to correct for them as they are
comparable in strength to flux densities of WINQSEs. As described in Sec. 3,
an additional step of forcing the median disc integrated solar flux density to
a constant value has been introduced in the calibration process to ensure a
constant fluxscale throughout the observing duration.

  2. 2.

During the calibration process, a varying number of antennas can get flagged
for different time slices even at the same frequency. This can lead to small
differences between the psfs for different images. To ensure that this does
not affect the flux density estimated, all of the images were smoothed to the
same coarser resolution for this work, as described in Sec. 3.

###  5.2 Comparison with M20

The results from this improved methodology confirm the key properties of
WINQSEs reported by M20 – weak flux densities, impulsive nature and ubiquitous
presence on the quiet Sun during this period of very low solar activity.
However, there are some significant quantitative differences which need some
discussion.

####  5.2.1 $\Delta F/F$ distribution of WINQSEs

This work highlights that the distribution of $\Delta F/F$, can change
significantly between observations. At the low $\Delta F/F$ end, the
distributions found by M20 tend to flatten out rather than drop, as seen in
Fig. 3. At the high $\Delta F/F$ end they show clear power law tails spanning
more than an order of magnitude which are clearly missing in this work.

For a given frequency, M20 plotted all $\Delta F/F$ values for all times and
all the psf sized regions when presenting these histograms, irrespective of
the signal-to-noise of the detection. Hence the low end of the $\Delta F/F$
distribution was dominated primarily by noise fluctuations. In contrast, we
use an image noise based threshold on $\Delta F$ to limit ourselves only to
reliable detections of WINQSEs. This filters out the low signal-to-noise
detections and explains the discrepancy between these works at the low end of
the $\Delta F/F$ histograms. To substantiate this, we plot in Fig. 7 the data
from an example frequency from M20 with the same thresholding technique as
used here. It is evident that the low $\Delta F/F$ end of the histogram is now
similar to the ones shown in Fig. 3, though it retains the powerlaw behaviour
shown in M20 at the large $\Delta F/F$ end. The powerlaw index estimated is
consistent with that reported in M20. Similar results were obtained at other
frequencies as well.

Though M20 had attempted to avoid contamination from the only active region
present on the Sun during their observations, we find that despite their
efforts the high $\Delta F/F$ tail observed in the histograms obtained in M20
arises due to this lone active region. To substantiate this, the left panel of
Fig. 8 shows the location of all WINQSEs detected in M20 at 132 MHz, overlaid
on an example radio image at the same frequency. The regions are colored blue
(red) if they host at least one WINQSE with $\Delta F/F>0$ ($>1$). It is
evident that while WINQSEs are distributed all over the Sun, those with
$\Delta F/F>1$ are clustered around the active region and not a single one of
them is present far away from it.

Active regions have very complex magnetic fields and are sites of continuous
magnetic reconnections (e.g. Mondal & Oberoi, 2021). The complexity of the
magnetic field is expected to decrease smoothly as the distance from the
active region increases. Consequently, the strength and the number of magnetic
reconnections taking place would also decrease in a similar smooth manner. We
hypothesise this to be the reason for regions showing the presence of strong
WINQSEs in M20 to lie close to the active region. Additionally, such a region
is expected not only to produce strong WINQSEs, but also to give rise to a
larger number of weaker WINQSEs, as compared to regions far away from the
active region. Hence we hypothesise this to be the reason behind the presence
(absence) of powerlaw behaviour in the $\Delta F/F$ distribution of the
detected WINQSEs in M20 (this work). To corroborate this hypothesis, the right
panel of Fig. 8 presents the $\Delta F/F$ distribution after removing the
contribution of the regions shown in red in the left panel. The red line shows
the straight line fitted to the same range of $\Delta F/F$ used for the fit in
Fig. 7. It is evident that the powerlaw is no longer a good representation of
the distribution. It is also observed that there is hump-like feature around
$\Delta F/F\approx 0.5$. The signature of this hump is also observed at a
similar location in Fig. 7 in this work and Fig 4 of M20. Understanding the
exact details of this feature is beyond the scope of the present work. This
observation of the transition from a powerlaw behaviour towards a more
lognormal-like nature suggests that as the solar activity increases, more
active regions are observed, leading to a larger number of stronger WINQSEs.
Hence with increasing solar activity a transition from a lognormal nature to a
powerlaw behaviour in the $\Delta F/F$ distribution can be expected. This
prediction should be tested in future with data spanning a large part of a
solar cycle.

Note that the peaks in the $\Delta F/F$ histograms lie at lower values in the
current data as compared to M20. We attribute this to a combination of
exceptionally quiet solar conditions and lower image noise characteristics
delivered by the improved calibration employed here.

Figure 7: The data used here correspond to the 132 MHz data presented in M20.
The $\Delta F/F$ histogram using data from M20 following the same procedure as
that followed here. The red line shows the powerlaw fit to the data beyond the
peak. The powerlaw index is mentioned in the figure.

Figure 8: Left panel: The blue filled circles show the locations of all the
psf sized regions used by M20 with WINQSEs detections. The regions which were
determined by M20 to be affected by active region emission and excluded from
their analysis have been excluded here as well. The red filled circles show
the locations where the $\Delta F/F$ of the detected WINQSEs exceeds 1. The
contours correspond to 0.01, 0.02, 0.04, 0.08, 0.16, 0.32, 0.64 times the peak
of an example image at 132 MHz. Histogram made using the same data as
presented in Fig. 7 after removing contribution of regions which show at least
one $\Delta F/F>1$.

####  5.2.2 Fractional occupancy of WINQSEs

The median fractional occupancies obtained from the data presented here are
0.017, 0.021, 0.024 and 0.045 at 120, 128, 136 and 146 MHz respectively. In
comparison, at similar frequencies of 120 and 132 MHz, M20 had found higher
median fractional occupancies of 0.07 and 0.06 respectively. There are at
least four known factors which contribute to a higher fractional occupancy for
M20. The first is the improved analysis methodology used here which rigorously
rejects low signal-to-noise detections. This was not the case with M20 and
hence some spurious noise peaks can get counted as WINQSEs inflating the true
fractional occupancy number. Secondly, M20 computed a single median for the
entire duration, in contrast to the piece-wise median estimated here. The
appendix shows the results obtained using a single median and it is evident
that the results thus obtained are closer to that obtained by M20 (3-5%). The
third important reason is the more active state of the Sun during M20, some
evidence and arguments for contributions from which have already been
presented.

###  5.3 Comparison with S22

As mentioned earlier, S22 have recently reported the detection of WINQSEs
using a very different technique, namely the visibility subtraction technique.
For completeness, we briefly describe the technique first and then present a
comparison of their technique and the results with the present work.

Similar to the technique used here, S22 also assumes that there is a slowly
varying background solar radio emission on which are the superposed the
impulsive WINQSEs. While we separate the slowly varying and impulsive parts of
emission in the image domain by estimating and subtracting a median image, S22
does so in the visibility domain. For each baseline, they compute the median
complex visibility over a 15 second running window and assume that the median
represents the visibility from the slowly varying Sun for that baseline for
the instant centered on that time interval. In principle, subtraction in the
visibility domain and subtraction in the image domain should be identical due
to the Fourier relationship between them. However, in practice, there can be
differences due to non-linear nature of deconvolution which is used to produce
the radio images. The subtraction in visibility domain is expected to perform
better if there are significant deconvolution errors in the images. A
disadvantage of the visibility subtraction technique, however, is that the
Fourier component sampled by a given baseline changes as a function of time.
For a source with the complex morphology like the Sun, the emission morphology
can itself give rise to variations in the observed visbilities on a given
baseline even if the source were time invariant. To limit the magnitudes of
these intrinsic changes in measured visibilities, S22 used a rather short
duration of 15 s over which to compute the median, while our technique was
able to use the full 4 min duration of the scan to compute the median image.
An implication is that S22 would miss the longer duration events. While this
will bias any studies of the temporal width distribution done using this
technique, given the very steep slope of this distribution, however, the
number of events missed will be a very small fraction of the total number.

The solar images used in this work were of high quality and no deconvolution
errors or artefacts were evident. In addition to the excellent uv-coverage,
precise calibration performed using AIRCARS (Mondal et al., 2019) and the
absence of any strong source on the solar disc, all have contributed to
ensuring a high fidelity images. This has allowed us to use the image
subtraction technique with a high confidence and enabled the detection of
these WINQSEs.

###  5.4 Need for a better WINQSE detection technique

A shortcoming of both this work and M20 stems from the way WINQSEs are
detected. We have searched for WINQSEs in discrete psf-sized regions and
classified the transients as true detections if their flux density exceeds
some threshold. If a WINQSE is larger than the psf, then it will contribute
flux density to multiple regions and hence can lead to over-counting of the
number of WINQSEs. On the other hand, since the psf is circular/elliptical in
nature, the regions used here and in M20 do not form a closed packed
structure. Hence some part of the Sun is not covered by these regions and
could lead to some WINQSEs being missed. Additionally, when the spatial
location of a WINQSE does not coincide with the centre of some psf sized
region, the observed increase in flux density of the region due to the WINQSE
is lower than its true flux density. This will lead to instances when the
observed flux density will fall short of the detection threshold due to this
artificial reduction. Due to these effects, it is non-trivial to estimate the
true number of WINQSEs present in the data. To address this shortcoming, we
have recently developed a machine learning based technique to reliably
identify and characterize WINQSEs in solar radio images. This work will be
presented in an upcoming paper (Bawaji et al., 2022).

An additional shortcoming of both these works stems from the way the median is
computed. M20 had computed a single median map for the entire time duration
studied, whereas here we have computed a median map for every five minutes of
data. Both of these methods have their own strengths and weaknesses. For
example, while calculating an independent median for every five minutes of
data, can reduce some systematic calibration artefacts, if present, it can
also produce spurious WINQSEs at the five minute data boundaries if for that
data chunk and pixel, the lightcurve has a smoothly decreasing or increasing
trend. We have repeated the entire analysis using a single median map and, as
expected, we find that the number of WINQSEs found in this analysis is
significantly larger than when done using independent medians for each data
chunk. We are not aware of any systematic effects in the data which can change
the flux densities in different regions of the Sun while maintaining the
integrated flux at all times. Choosing to err on the side of caution, we have
presented the results corresponding to the medians computed every five minutes
in the preceding sections. The statistical properties of WINQSEs obtained from
the single median analysis are very similar and are presented in the Appendix.
While the results are consistent, we realize the need to improve WINQSEs
detection strategy beyond what has been employed here. Work on this front is
already underway.

###  5.5 A dynamic middle corona

Radio observations at the small range of frequencies used here probe coronal
height of $\sim 1.4\ (2.0)\ R_{\odot}$, under the standard assumption that the
emission arises due to plasma emission mechanism at the fundamental (harmonic)
and using well established coronal electron density models. The observations
in M20, Sharma et al. (2022) and this work imply that energetically weak
reconnections are ever-present even in the quiet Sun regions at these coronal
heights. These reconnections accelerate electrons which then emit in the radio
band. Due to instrumental limitations at high energies (EUV and X-rays), these
coronal heights are typically accessible only at metrewave radio frequencies.
The bulk of the radio studies in the past have focused on strong events, again
primarily due to the instrumental limitations in radio bands. So there have
been rather few investigations of weak or level emissions from the middle
coronal heights. As more capable instruments are becoming available, the
evidence for middle coronal heights to also be a dynamic place is slowly
building. Recently, gradual bulk flows originating in the middle corona have
been observed using the Solar Ultraviolet Imager (SUVI) onboard the
Geostationary Operational Environmental Satellite (Seaton et al., 2021). The
authors suggest that this can happen due to magnetic reconnection processes
higher in the corona. While a direct comparison between these EUV emissions
and WINQSEs is beyond the scope of this work, these disparate observations
paint a consistent picture of a highly dynamic middle solar corona even during
extremely quiet times.

##  6 Conclusion

M20 was the first work to present ubiquitous detection of WINQSEs. Due to the
significant potential implications of this work on coronal heating, and
presence of nonthermal particles in the solar corona in general, it is very
important to verify their presence under a variety of solar conditions, take a
closer look at the analysis and characterise them better. With this
motivation, this work forms the next step in this progression. We have
examined data from extremely quiet solar conditions, improved on the
methodology used in M20 and presented robust detection of large numbers of
WINQSEs. This work corroborates the key findings of M20 – their ubiquity and
impulsive nature. While significant differences between the distribution of
$\Delta F/F$ presented here and M20 are observed, we show that they can be
justifiably attributed to differences in methodology and variations in the
level of solar activity. Parallel to this effort, independent techniques have
been developed to identify and characterise WINQSEs (Bawaji et al., 2022;
Sharma et al., 2022) and they all paint a consistent picture. Together these
works place the presence of WINQSEs on a firm pedestal. We hope that this will
provide the motivation to the wider community to understand them in greater
detail.

This scientific work makes use of the Murchison Radio-astronomy Observatory,
operated by CSIRO. We acknowledge the Wajarri Yamatji people as the
traditional owners of the Observatory site. Support for the operation of the
MWA is provided by the Australian Government (NCRIS), under a contract to
Curtin University administered by Astronomy Australia Limited. We acknowledge
the Pawsey Supercomputing Centre which is supported by the Western Australian
and Australian Governments. We thank the anonymous referee for their comments,
which have improved this paper significantly. SM acknowledges partial support
by USA NSF grant AGS-1654382 to the New Jersey Institute of Technology. AB and
DO acknowledge support of the Department of Atomic Energy, Government of
India, under the project no. 12-R&D-TFR-5.02-0700. As discussed in Sec. 5.4,
the number of WINQSEs detected depends on the choice of the median map used.
Results presented in the paper are based on the use of an independent median
computed for each five minute dataset. For completeness and comparison, here
we present results (analogs of Figs. 3–6) obtained using a single median
computed over the entire observing duration as was done by Mondal et al.
(2020). As expected the number of WINQSEs detected here are larger than using
a median computed every five minutes.

Figure 9: Histogram of $\Delta F/F$ using a median computed over the entire
duration of the data is shown using the blue circles. The red line shows the
lognormal fit to these data. The observation frequencies are written in each
panel. Figure 10: Temporal width distribution of the detected WINQSEs using a
median computed over the entire duration of the data. Figure 11: Wait time
distribution of the detected WINQSEs using a median computed over the entire
duration of the data. The red and black lines show the wait time till which
the uncertainties due to incomplete sampling is less than 10% and 20%
respectively. Figure 12: Fractional occupancy distribution of the detected
WINQSEs using a median computed over the entire duration of the data.

## References

  * Aschwanden & McTiernan (2010) Aschwanden, M. J., & McTiernan, J. M. 2010, ApJ, 717, 683, doi: 10.1088/0004-637X/717/2/683
  * 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
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Bain & Fletcher (2009) Bain, H. M., & Fletcher, L. 2009, A&A, 508, 1443, doi: 10.1051/0004-6361/200911876
  * Bawaji et al. (2022) Bawaji, S., Alam, U., Mondal, S., & Oberoi, D. 2022, in Astronomical Society of the Pacific Conference Series, Vol. 532, Astronomical Society of the Pacific Conference Series, ed. J. E. Ruiz, F. Pierfedereci, & P. Teuben, 211.  https://arxiv.org/abs/2103.05371
  * Che (2018) Che, H. 2018, in Journal of Physics Conference Series, Vol. 1100, Journal of Physics Conference Series, 012005, doi: 10.1088/1742-6596/1100/1/012005
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Kansabanik et al. (2022a) Kansabanik, D., Mondal, S., Oberoi, D., Biswas, A., & Bhunia, S. 2022a, ApJ, 927, 17, doi: 10.3847/1538-4357/ac4bba
  * Kansabanik et al. (2022b) Kansabanik, D., Oberoi, D., & Mondal, S. 2022b, ApJ, 932, 110, doi: 10.3847/1538-4357/ac6758
  * Kivelä & Porter (2015) Kivelä, M., & Porter, M. A. 2015, Phys. Rev. E, 92, 052813, doi: 10.1103/PhysRevE.92.052813
  * Kontar et al. (2017) Kontar, E. P., Yu, S., Kuznetsov, A. A., et al. 2017, Nature Communications, 8, 1515, doi: 10.1038/s41467-017-01307-8
  * Krucker et al. (2011) Krucker, S., Kontar, E. P., Christe, S., Glesener, L., & Lin, R. P. 2011, ApJ, 742, 82, doi: 10.1088/0004-637X/742/2/82
  * Lonsdale et al. (2009) Lonsdale, C. J., Cappallo, R. J., Morales, M. F., et al. 2009, Proceedings of the IEEE, 97, 1497, doi: 10.1109/JPROC.2009.2017564
  * McLean & Labrum (1985) McLean, D. J., & Labrum, N. R. 1985, Solar radiophysics : studies of emission from the sun at metre wavelengths
  * McMullin et al. (2007) McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in Astronomical Society of the Pacific Conference Series, Vol. 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell, 127
  * Melrose (1980) Melrose, D. B. 1980, Space Sci. Rev., 26, 3, doi: 10.1007/BF00212597
  * Mohan (2021) Mohan, A. 2021, A&A, 655, A77, doi: 10.1051/0004-6361/202142029
  * Mohan et al. (2019) Mohan, A., Mondal, S., Oberoi, D., & Lonsdale, C. J. 2019, ApJ, 875, 98, doi: 10.3847/1538-4357/ab0ae5
  * Mondal (2021) Mondal, S. 2021, Sol. Phys., 296, 131, doi: 10.1007/s11207-021-01877-3
  * Mondal et al. (2019) Mondal, S., Mohan, A., Oberoi, D., et al. 2019, The Astrophysical Journal, 875, 97, doi: 10.3847/1538-4357/ab0a01
  * Mondal & Oberoi (2021) Mondal, S., & Oberoi, D. 2021, ApJ, 920, 11, doi: 10.3847/1538-4357/ac1076
  * Mondal et al. (2020) Mondal, S., Oberoi, D., & Mohan, A. 2020, The Astrophysical Journal, 895, L39, doi: 10.3847/2041-8213/ab8817
  * Parker (1988) Parker, E. N. 1988, ApJ, 330, 474, doi: 10.1086/166485
  * Pauluhn & Solanki (2007) Pauluhn, A., & Solanki, S. K. 2007, Astronomy and Astrophysics, 462, 311, doi: 10.1051/0004-6361:20065152
  * Pick & Vilmer (2008) Pick, M., & Vilmer, N. 2008, A&A Rev., 16, 1, doi: 10.1007/s00159-008-0013-x
  * Reid (2020) Reid, H. A. S. 2020, Frontiers in Astronomy and Space Sciences, 7, 56, doi: 10.3389/fspas.2020.00056
  * Reid & Kontar (2017) Reid, H. A. S., & Kontar, E. P. 2017, A&A, 606, A141, doi: 10.1051/0004-6361/201730701
  * Seaton et al. (2021) Seaton, D. B., Hughes, J. M., Tadikonda, S. K., et al. 2021, Nature Astronomy, 5, 1029, doi: 10.1038/s41550-021-01427-8
  * Sharma et al. (2022) Sharma, R., Oberoi, D., Battaglia, M., & Krucker, S. 2022, ApJ, 937, 99, doi: 10.3847/1538-4357/ac87fc
  * Sturrock (1964) Sturrock, P. A. 1964, Type III Solar Radio Bursts, Vol. 50, 357
  * Tajfirouze & Safari (2012) Tajfirouze, E., & Safari, H. 2012, ApJ, 744, 113, doi: 10.1088/0004-637X/744/2/113
  * Takakura (1967) Takakura, T. 1967, Sol. Phys., 1, 304, doi: 10.1007/BF00151359
  * Tingay et al. (2013) Tingay, S. J., Goeke, R., Bowman, J. D., et al. 2013, Publications of the Astronomical Society of Australia, 30, doi: 10.1017/pasa.2012.007
  * Upendran & Tripathi (2021) Upendran, V., & Tripathi, D. 2021, ApJ, 916, 59, doi: 10.3847/1538-4357/abf65a
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Wayth et al. (2018) Wayth, R. B., Tingay, S. J., Trott, C. M., et al. 2018, PASA, 35, e033, doi: 10.1017/pasa.2018.37
  * Wild & Smerd (1972) Wild, J. P., & Smerd, S. F. 1972, ARA&A, 10, 159, doi: 10.1146/annurev.aa.10.090172.001111
  * Wild et al. (1963) Wild, J. P., Smerd, S. F., & Weiss, A. A. 1963, ARA&A, 1, 291, doi: 10.1146/annurev.aa.01.090163.001451


compat=1.0.0

# Singlet Scalars as Dark Matter and the Muon $(g-2)$ Anomaly

Bastián Díaz Sáez<EMAIL_ADDRESS>Physik-Department, Technische
Universität München, James-Franck-Straße, 85748 Garching, Germany  Karim
Ghorbani<EMAIL_ADDRESS>Physics Department, Faculty of Sciences,
Arak University, Arak 38156-8-8349, Iran

###### Abstract

We explore a simple model containing two singlet scalars as dark matter
candidates and a vector-like lepton. Evading collider and dark matter
constraints, the model is able to accommodate the correct dark matter relic
abundance and the $(g-2)_{\mu}$ anomaly when the masses of the three new
fields are around 100 GeV and the values of the minimum amount of required
couplings are order one.

##  1 Introduction

The reported combined measurement of the anomalous magnetic moment of the muon
$a_{\mu}\equiv\frac{g-2}{2}$ by BNL E821 [1] and FNAL E989 [2] experiments
yields a 4.2$\sigma$ discrepancy from the Standard Model (SM) prediction [3,
4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24]

$\displaystyle\Delta a_{\mu}=a_{\mu}^{\text{exp}}-a_{\mu}^{\text{SM}}=(251\pm
59)\times 10^{-11}.$ (1)

This hint of new physics beyond the Standard Model (BSM) has motivated
different extensions to the SM to account for it. In particular, new
interactions involving the leptonic sector along with dark matter (DM)
candidates seem an appealing and exiting avenue [25, 26, 27, 28, 29, 30, 31,
32, 33, 34, 35, 36, 37, 38, 39, 40, 41]. One of the simplest extensions to the
SM to account for DM and the $(g-2)_{\mu}$ anomaly involves the introduction
of two new fields, i.e. one real scalar plus one vector-like lepton (VLL),
both odd under the same discrete symmetry $Z_{2}$ ($S_{1}\rightarrow-S_{1}$
and $\psi\rightarrow-\psi$). However, this set-up has been shown to give an
under-abundant relic density [33]. Extensions with three fields have also been
explored [28, 33], and in particular, it has been shown that extending the
minimal scenario (singlet scalar plus a VLL) by including a $SU(2)_{L}$ scalar
doublet may accommodate $(g-2)_{\mu}$ and DM constraints, provided that the
scalar DM masses be above 500-600 GeV, and in the TeV scale when the mixed
scalars are degenerated [33].

In this work we explore a simpler scenario in which we extend the minimal
model (i.e. one singlet scalar $S_{1}$ and one VLL) by adding another singlet
scalar $S_{2}$, which along with $S_{1}$ both play the role of DM. The model
is able to explain dark matter and $(g-2)_{\mu}$ at the EW scale.
Multicomponent singlet scalar DM has been an interesting subject to explore
[42, 43, 44], and in particular it has been shown to account for other event
excess, such as the cosmic-ray positron-to-electron ratio [45].

The paper is organized in the following way. In section 2 we introduce the
model and its freeze-out mechanism. In section 3 we present the most relevant
constraints on the model. In section 4 we present different scenarios to
account for $(g-2)_{\mu}$ and DM, and finally in section 5 a discussion with
the conclusions.

##  2 Set-up

###  2.1 Lagrangian

We introduce two real singlet scalars $S_{1}$ and $S_{2}$, and a fermion
$\psi$ being singlet under $SU(2)_{L}$ with $Y=-1$. Two discrete symmetries
are imposed $Z_{2}\times Z_{2}^{\prime}$, where the new fields transform as

$\displaystyle Z_{2}$ $\displaystyle:$ $\displaystyle
S_{1},S_{2},\psi\rightarrow-S_{1},S_{2},-\psi,$ (2) $\displaystyle
Z_{2}^{\prime}$ $\displaystyle:$ $\displaystyle S_{1},S_{2},\psi\rightarrow
S_{1},-S_{2},\psi.$ (3)

In this way, the Lagrangian is given by

$\displaystyle\mathcal{L}$ $\displaystyle\supset$ $\displaystyle
m_{\psi}\bar{\psi}\psi+(\kappa_{1}\bar{\mu}_{R}S_{1}\psi_{L}+h.c.)$
$\displaystyle+$
$\displaystyle\sum_{i=1,2}\left(\frac{\mu_{i}^{2}}{2}S_{i}^{2}+\frac{\lambda_{Hi}}{2}|H|^{2}S_{i}^{2}+\frac{\lambda_{i}}{4!}S_{i}^{4}\right)+\frac{\lambda_{12}}{4}S_{1}^{2}S_{2}^{2}$

All the parameters in the model are real. For the analysis of this paper the
quartic couplings, $\lambda_{i}$, play no role, so we can set $\lambda_{i}$ as
arbitrary values in our computations. Additionally, none of the scalar fields
get vev, and due to the two discrete symmetries, both $S_{1}$ and $S_{2}$ are
stable. After electroweak symmetry breaking (EWSB), the masses of both singlet
scalars at tree level are given by
$m_{i}^{2}=\mu_{i}^{2}+\lambda_{Hi}v_{H}^{2}/2$, with $i=1,2$.

###  2.2 Relic abundance

As the model contains two DM components, the total relic abundance is given by
the sum of its partial contributions, i.e. $\Omega
h^{2}=\Omega_{1}h^{2}+\Omega_{2}h^{2}$, whose value must be contrasted with
the Planck measurement $\Omega_{c}h^{2}=0.12$ [46]. In order to achieve the
correct relic abundance and at the same time respecting the new measurements
on $(g-2)_{\mu}$, in our framework it requires to have the following mass
hierarchy

$\displaystyle m_{2}\lesssim m_{1}<m_{\psi}.$ (5)

The first inequality is necessary in such a way that DM conversions of the
type $S_{1}S_{1}\rightarrow S_{2}S_{2}$ be effective in the thermal decoupling
history, yielding $\Omega_{2}>\Omega_{1}$, and in this way $\Omega_{2}$
covering the missing DM in the minimal scenario of two-fields scalar plus VLL
[33]. The second inequality is necessary in order to have $S_{1}$ stable,
otherwise it could decay into the pair $\psi^{\pm}\mu^{\mp}$.

Introducing $x=\mu/T$, with $\mu=m_{1}m_{2}/(m_{1}+m_{2})$ and $T$ the SM
temperature, and $Y_{i}=n_{i}/s$, the coupled Boltzmann equations for the
system of the two singlet scalars in thermal equilibrium in the early universe
are given by

$\displaystyle\begin{split}\dfrac{dY_{1}}{dx}=&-\lambda_{11XX}\cdot[Y_{1}^{2}-Y_{1e}^{2}]-\lambda_{1122}\cdot[Y_{1}^{2}-r_{12}^{2}Y_{2}^{2}]\\\
\dfrac{dY_{2}}{dx}=&-\lambda_{22XX}\cdot[Y_{2}^{2}-Y_{2e}^{2}]+\lambda_{1122}\cdot[Y_{1}^{2}-r_{12}^{2}Y_{2}^{2}],\end{split}$
(6)

with $r_{ij}(x):=Y_{ie}(x)/Y_{je}(x)$ and

$\displaystyle\lambda_{iijj}(x):=\dfrac{\langle\sigma_{iijj}v\rangle(x)\cdot
s(T)}{x\cdot H(T)},\qquad\text{for}\leavevmode\nobreak\ i,j=1,2,X.$ (7)

where $\langle\sigma_{iijj}v\rangle(x)$ is the Maxwell-Boltzmann thermal
average annihilation cross section times the relative velocity, $s(T)$ and
$H(T)$ the entropy density and Hubble rate, respectively, and $X$ stands for
the SM particles, including the Higgs boson (for more details about the
Boltzmann equations and its components see [47]). One important point here is
that the Boltzmann equations in eq. (6) do not include coannihilation
processes, e.g. $S_{1}+\psi^{\pm}\rightarrow\mu^{\pm}+V$, with $V=\gamma,Z$,
due to the fact that they become effective provided that
$m_{\psi}^{2}/m_{1}^{2}\lesssim 1.2$ [48]. We will see in the next section
that this mass region is excluded by soft lepton searches, and therefore
making coannihilation processes ineffective in the parameter space of our
interest.

To illustrate the mechanism of DM production (also known as assisted-freeze-
out mechanism [49, 47]), in Fig. 1 we exemplify the thermal evolution of the
yields of the two scalar singlets. We choose the three masses and $\kappa_{1}$
in the typical ball park to account for $\Delta a_{\mu}$: $m_{\psi}=150$ GeV,
$m_{1}=130$ GeV, $m_{2}=125$ GeV and $\kappa_{1}=2.8$. For simplicity we set
$\lambda_{H1}=\lambda_{H2}=0$. The yields are shown for $\lambda_{12}=0.1$
(dashed lines) and $\lambda_{12}=1$ (solid lines). As it can be seen in the
plot, after $S_{1}$ (solid blue lines) decouples from the thermal bath (dotted
lines), it continues annihilating into $S_{2}$ (red lines), also known as DM
conversions, reducing its abundance even more. As a consequence, the scenario
always predicts $\Omega_{1}\leavevmode\nobreak\ \ll\leavevmode\nobreak\
\Omega_{2}$. As a matter of comparison, we show the thermal yield of a single
DM scalar in a VLL portal (see Sec. 1) with the green dashed line, making
explicit its sub-abundance [33].

Figure 1: Relic yields evolution of $S_{1}$ (blue curves) and $S_{2}$ (red
curves) as a function of $x$. The parameters are $(m_{1},m_{2},m_{\psi})$ =
(130, 125, 150) [GeV], $\lambda_{H1}=\lambda_{H2}=0$ and $\kappa_{1}=2.8$. The
dashed lines consider $\lambda_{12}=0.1$, whereas for the solid ones
$\lambda_{12}=1$. The green dashed line is the yield of single DM component
$S_{1}$ in a VLL portal, as described in Sec. 1.

One important aspect of this type of mechanism is the relevance of the mass
shift $\Delta m\equiv m_{1}-m_{2}$, which must not be bigger than a few GeV
[47]. As it can be noted from the second line of eq. (6), the relevant term
depending on the mass shift is

$\displaystyle\frac{dY_{2}}{dx}\supset-\lambda_{1122}r_{12}^{2}Y_{2}^{2}\sim-\lambda_{12}^{2}e^{-x\frac{\Delta
m}{\mu}}Y_{2}^{2}.$ (8)

This implies that the decreasing rate of $Y_{2}$ is not only regulated by
$\lambda_{12}$, but it also depends strongly on $\Delta m$ and $\mu$. This is,
if $\Delta m\rightarrow 0$, there is no exponential suppression, then
enhancing the decreasing of $Y_{2}$, and on the other hand, higher scalar
masses (bigger $\mu$) reduces $Y_{2}$ more efficiently. This last aspect will
be part of the analysis in Sec. 4.

Finally, we solved the system of Boltzmann equations 6 with our own code in
Python, and cross-checked our results using Micromegas code [50].

##  3 Constraints

On the scenarios in the present work we take into account the following
constraints:

  * •

Dark Matter. First, we set that the total relic predicted by the model $\Omega
h^{2}$ must fulfill the Planck measurement $\Omega_{c}h^{2}=0.12$ [46], with a
tolerance of $\pm$10% accounting for uncertainties in our computation.
Secondly, direct detection (DD) bounds for DM candidates with masses of
$\mathcal{O}(100)$ GeV are set by Xenon1T [51]. Two types of portals may be
present in the model, and each signal must be weighted by the relative
abundance of the corresponding DM component. DM-nucleon interactions via the
VLL portal start at two-loops, but this contribution turns out to be well
below Xenon1T bounds [52, 53]. In the case of Higgs portal interactions the
spin-independent DM-nucleon scattering amplitude is given by [54, 55]

$\displaystyle\sigma_{\text{SI}}=\frac{\lambda_{Hi}^{2}f_{N}}{4\pi}\frac{\mu^{2}m_{n}^{2}}{m_{h}^{4}m_{i}^{2}},\quad
i=1,2.,$ (9)

where $\mu=m_{n}m_{i}/(m_{n}+m_{i})$, $f_{N}=0.3$ and $m_{n}$ the nucleon
mass. Finally, the leading indirect detection (ID) signals today are
$S_{1}S_{1}\rightarrow\mu\mu(+\gamma)$ or $S_{1}S_{1}(S_{2}S_{2})\rightarrow
XX$, with $X=Z,W^{\pm},b(\bar{b})$ from the VLL and Higgs portal,
respectively. We used the most recent Fermi-Lat bounds given in [56].

  * •

Collider. LEP/LHC constrains the mass of the charged VLL to be above $\sim$100
GeV [57]. Furthermore, the Higgs decay into a pair of muon/antimuon may give
rise to important contributions at one-loop level. Considering the measured
cross section for the SM Higgs boson production, the upper limit at 95% C.L.
is given by $\text{Br}(h\rightarrow\mu^{+}\mu^{-})\leavevmode\nobreak\
<4.7\cdot\leavevmode\nobreak\ 10^{-4}$ [58, 59]. Taking into account only
$\lambda_{H1}$ coupling (equivalently for $\lambda_{H2}$), the one-loop decay
width is given by

$\displaystyle\Gamma(h\rightarrow\mu^{+}\mu^{-})=\frac{\kappa_{1}^{4}\lambda_{H1}^{2}v_{H}^{2}m_{h}}{4\pi
m_{\psi}^{2}}\frac{m_{\mu}^{2}}{m_{\psi}^{2}}|F(x,y)|^{2},$ (10)

with $F(x,y)=\frac{y}{x(y-1)}\log y+\frac{y-x}{x(x+y+1)}\log(y-x)$, in which
$x=m_{h}^{2}/m_{\psi}^{2}$ and $y\equiv m_{1}^{2}/m_{\psi}^{2}$. Finally,
searches for SUSY long-lived particles and also the CMS searches for the
compressed spectra of soft leptons, set limits on the combination
$(m_{\psi},m_{1})$. Here we take the results of [33], where the bounds found
there apply directly to our scenario.

  * •

Flavor Physics. We take the result in 1 as the reference deviation in the muon
anomalous magnetic moment. The leading contribution to $(g-2)_{\mu}$ in this
model comes from the one-loop diagram depicted in Fig. 2, and it is given by
[60]

$\displaystyle\Delta
a_{\mu}=\frac{\kappa_{1}^{2}}{96\pi^{2}}\frac{m_{\mu}^{2}}{m_{\psi}^{2}}f\left(\frac{m_{S_{1}}^{2}}{m_{\psi}^{2}}\right),$
(11)

where $f(t)=(2t^{3}+3t^{2}-6t^{2}\log t-6t+1)/(t-1)^{4}$, which is positive
for any $t$, and $f(0)=1$.

$\gamma$$\psi$$\bar{\psi}$$\mu$$\bar{\mu}$$S_{1}$ Figure 2: One-loop
contributions to $(g-2)_{\mu}$ in this work.

##  4 Analysis

In the following we show how the model may account for the correct relic
abundance and $(g-2)_{\mu}$, evading the constraints previously described.
Although the parameter space of the model is

$\displaystyle\left(m_{\psi},m_{1},m_{2},\kappa_{1},\lambda_{12},\lambda_{H1},\lambda_{H2}\right)$
(12)

we start the analysis neglecting the Higgs portal, and afterwards we discuss
the effects of it when the portal is not neglected.

Figure 3: Results in the scenario with $\lambda_{H1}=\lambda_{H2}=0$, scanning
over the masses of the new fermion and the heaviest scalar $S_{1}$, for
$\kappa_{1}=2.5$ (plot on the top) and $\kappa_{1}=3.5$ (plot on the bottom).
The region that can explain $\Delta a_{\mu}$ to within 1$\sigma$ is colored
green. The blue and red curves indicate those points that produce the correct
relic abundance, with the scalar mass difference $\Delta m$ and $\lambda_{12}$
indicated in the legend. The blue region is not allowed by instability of
$m_{1}$, and the yellow and grey regions are LHC 13 TeV and soft leptons
(s.l.) searches obtained in [33] (for details see the experimental references
therein).

###  4.1 Higgs portals closed

We start assuming that the Higgs portals are sub-dominant, i.e.
$(\lambda_{H1},\lambda_{H2})\ll 1$. In this case, the model becomes very
simple since only the first four parameters of 12 are relevant, and the
constraints presented in the last section reduce to a few ones: relic
abundance, $(g-2)_{\mu}$ and some collider constraints. Then, for simplicity
we assume that $\lambda_{H1}=\lambda_{H2}=0$.

In Fig. 3 we show some predictions of the model projected on the plane
$(m_{\psi},m_{1})$ in which all the constraints are fulfilled and
$(g-2)_{\mu}$ is accounted, for $\kappa_{1}=2.5$ (plot on the top) and
$\kappa_{1}=3.5$ (plot on the bottom). The white region is still allowed by
the existing constraints, whereas the yellow and grey regions are already
ruled out by the collider bounds. The blue region is not allowed in this model
because it would make $S_{1}$ unstable
($S_{1}\rightarrow\psi^{\pm}+\mu^{\mp}$), and the purple one is where
coannihilations become effective. The green band accounts for $(g-2)_{\mu}$ at
1$\sigma$, with the black line near the middle corresponds to the central
value.

As it can be noted from the plot in the upper panel of Fig. 3, $m_{\psi}$ and
$m_{1}$ spanning several dozens of GeV with small $\Delta m$ and
$\lambda_{12}\approx 1$ enter the region allowed by the experimental
constraints fulfilling $(g-2)_{\mu}$. In particular, red curves with $\Delta
m=2$ GeV enter the green zone, with tiny deviations as $\lambda_{12}$ changes,
and as $\Delta m$ increases the model predictions deviate from the allowed
region, as it can be noted for $\Delta m=3$ GeV (blue curves). The strong
sensitivity of the model predictions to the growing of $\Delta m$ can be
understood due to the fact that in order to keep the correct relic abundance
at higher $\Delta m$, it is necessary to increment the reduced mass $\mu$ to
compensate the exponential suppression given in eq. (8). The sensitivity of
the model predictions to $\lambda_{12}$ peaks at high masses where the
exponential suppression in eq. (8) is less strong, therefore making the
abundance $Y_{2}$ being sensitive mainly to $\lambda_{12}^{2}$. Equivalently,
the plot in the lower panel of Fig. 3 considering $\kappa_{1}=3.5$ shows that
the viable mass shift picks up a little bigger value, as it is illustrated for
$\Delta m=5$ GeV (blue curves), with an increasing range for the masses
$m_{\psi}$ and $m_{1}$. Note that as a big portion of the coannihilation
region is excluded by LHC bounds (yellow region) and it practically does not
match the $(g-2)_{\mu}$ region (green one) we did not show the predictions of
the model in that region.

In this simple scenario, as the Higgs portal is closed, the leading
contributions to ID are
$S_{1}S_{1}\rightarrow\mu^{+}\mu^{-},\mu^{+}\mu^{-}\gamma,\gamma\gamma$.
However, in our scenarios $S_{1}$ is always a sub-dominant DM component,
typically $\Omega_{1}h^{2}\sim 10^{-3}$, then the resulting re-scaled average
annihilation cross section at zero velocity is suppressed by a factor of order
$\sim 10^{-6}$. We have checked that the resulting signals are some order of
magnitude below the upper bounds from Fermi-LAT.

###  4.2 Higgs portals opened

Now we consider the case that either one or two Higgs portals are not
negligible. In this way, here we sketch what deviations occur once the Higgs
portals are not negligible with respect to the previous scenario. Let us start
assuming that $\lambda_{H1}=\lambda_{H2}\neq 0$. It is not possible to
arbitrarily increase the values of the Higgs portal, because they not only
will deviate the relic abundance found in the previous scenario, but they
increase the DD signal. Without loss of generality, let us consider a point of
the parameter space that fulfills $(g-2)_{\mu}$ and LEP/LHC constraints:
$(m_{\psi},m_{1},m_{2})=(150,130,125)$ GeV, $\kappa_{1}=2.8$ and
$\lambda_{12}=1$. In Fig. 4(upper panel) we show the deviation in the relic
abundances as the Higgs portal couplings change according to
$\lambda_{H1}=\lambda_{H2}$.

Figure 4: Relic abundances as a function of Higgs portal parameters, with
$(m_{\psi},m_{1},m_{2})=(150,130,125)$ GeV, $\kappa_{1}=2.8$ and
$\lambda_{12}=1$. In the plot above we consider $\lambda_{H1}=\lambda_{H2}$,
whereas in below $\lambda_{H2}=0$. The blue and red curves are allowed by DD
bounds, whereas the grey one is the already ruled out by Xenon1T. The green
region is forbidden by Br$(h\rightarrow\mu^{+}\mu^{-})$ measurements, and the
dotted horizontal line is the measured relic abundance $\Omega_{c}h^{2}$.

For $\lambda_{H2}\gtrsim 0.01$, modification on $\Omega_{2}$ starts appearing
with respect to the case $\lambda_{H1}=\lambda_{H2}=0$ and almost
simultaneously the model enters the region in conflict with the DD bounds
(grey curve). Therefore, provided that the Higgs portal couplings take values
lower than $\mathcal{O}(0.01)$, no significant deviations with respect to what
we have found in Sec. 4.1 will appear in this case.

Since the scalar $S_{2}$ gives most of the DM relic abundance (i.e. in this
scenario we always have $\Omega_{1}\ll\Omega_{2}$) and therefore most sensible
to the DD bounds via the Higgs portal, another interesting option is to set
$\lambda_{H1}\neq 0$ keeping $\lambda_{H2}$ very small. In this case, higher
values for $\lambda_{H1}$ can be reached before the DD bounds start to be
relevant. In particular, for the same parameter space point described in the
previous paragraph, in Fig. 4(lower plot) $\lambda_{H1}\lesssim 0.1$ still
fulfills the relic abundance while evading DD bounds. Therefore, allowing the
Higgs portals coupling take non-vanishing values, DD bounds kill all the
significant deviations to the predictions obtained in Sec. 4.1 (Higgs portal
closed).

As a matter of complement, note that the scenario $\lambda_{H1}\neq 0$ and
$\lambda_{H2}=0$ is not unstable under radiative corrections due to the
following. At one-loop, the effective $S_{2}S_{2}h$ coupling in the limit in
which $\lambda_{H2}\ll 1$ is given by [47]

$\displaystyle\lambda_{H2}^{\text{eff}}=\lambda_{H2}+\frac{\lambda_{H1}\lambda_{12}}{16\pi^{2}}.$
(13)

As we have seen, the maximum value of $\lambda_{H1}$ allowed by the DD bounds
is $\mathcal{O}(0.1)$, then implying that at one-loop
$\lambda_{H2}^{\text{eff}}\approx 6.3\times 10^{-4}$, for
$\lambda_{12}\sim\mathcal{O}(1)$. As we have seen in this section, such small
values for $\lambda_{H2}^{\text{eff}}$ produce no deviation on the relic
abundance and evade easily the DD bounds.

##  5 Discussion and Conclusions

In this work we have shown that a simple unexplored extension to the SM
containing two real singlet scalars and one vector-like lepton not only
accounts for the measured DM relic abundance evading stringent experimental
constraints, but it can also fulfill the recent $4.2\sigma$ deviation in
$(g-2)_{\mu}$, to be tested in the forthcoming years with more data. Provided
that the mass of the new three states be around the Fermi scale and the mass
shift between the singlet scalars is not more than a few GeV, the
$(g-2)_{\mu}$ anomaly is perfectly achieved, with the minimal requirement that
the new couplings be of order one. We have seen that this set-up implies that
the inert scalar $S_{2}$ captures most of the relic abundance, whereas the
sub-leading component $S_{1}$ plays the role of enhancing the anomalous
magnetic moment of the muon at one-loop. Furthermore, we have seen that the
Higgs portal couplings play a secondary role in this framework, although they
can still take relatively sizable values in agreement with all the
experimental constraints. Finally, an additional comment about the possibility
that also $S_{2}$ couples to the muon via the VLL such that both scalars
participate in the $(g-2)_{\mu}$ (this can be done assuming
$\psi\rightarrow-\psi$ under both $Z_{2}$ and $Z_{2}^{\prime}$). In that case,
one can have $\kappa_{1}$ and $\kappa_{2}$ with slightly lower values than the
obtained ones in this work. That scenario is dangerous because the scalar with
the highest DM relic abundance (i.e. at least half of $\Omega_{c}h^{2}$) will
present strong indirect signals (e.g. bremsstrahlung), then entering in
conflict with the present ID bounds.

##  6 Acknowledgment

B.D.S has been founded by ANID (ex CONICYT) Grant No. 74200120.

## References

  * [1] Muon g-2 Collaboration, G. W. Bennett et al., “Final Report of the Muon E821 Anomalous Magnetic Moment Measurement at BNL,” Phys. Rev. D 73 (2006) 072003, arXiv:hep-ex/0602035.
  * [2] Muon g-2 Collaboration, B. Abi et al., “Measurement of the Positive Muon Anomalous Magnetic Moment to 0.46 ppm,” Phys. Rev. Lett. 126 no. 14, (2021) 141801, arXiv:2104.03281 [hep-ex].
  * [3] T. Aoyama, M. Hayakawa, T. Kinoshita, and M. Nio, “Complete Tenth-Order QED Contribution to the Muon $g-2$,” Phys. Rev. Lett. 109 (2012) 111808, arXiv:1205.5370 [hep-ph].
  * [4] T. Aoyama, T. Kinoshita, and M. Nio, “Theory of the Anomalous Magnetic Moment of the Electron,” Atoms 7 no. 1, (2019) 28.
  * [5] A. Czarnecki, W. J. Marciano, and A. Vainshtein, “Refinements in electroweak contributions to the muon anomalous magnetic moment,” Phys. Rev. D67 (2003) 073006, arXiv:hep-ph/0212229 [hep-ph]. [Erratum: Phys. Rev. D73, 119901 (2006)].
  * [6] C. Gnendiger, D. Stöckinger, and H. Stöckinger-Kim, “The electroweak contributions to $(g-2)_{\mu}$ after the Higgs boson mass measurement,” Phys. Rev. D88 (2013) 053005, arXiv:1306.5546 [hep-ph].
  * [7] M. Davier, A. Hoecker, B. Malaescu, and Z. Zhang, “Reevaluation of the hadronic vacuum polarisation contributions to the Standard Model predictions of the muon $g-2$ and ${\alpha(m_{Z}^{2})}$ using newest hadronic cross-section data,” Eur. Phys. J. C 77 no. 12, (2017) 827, arXiv:1706.09436 [hep-ph].
  * [8] A. Keshavarzi, D. Nomura, and T. Teubner, “Muon $g-2$ and $\alpha(M_{Z}^{2})$: a new data-based analysis,” Phys. Rev. D97 no. 11, (2018) 114025, arXiv:1802.02995 [hep-ph].
  * [9] G. Colangelo, M. Hoferichter, and P. Stoffer, “Two-pion contribution to hadronic vacuum polarization,” JHEP 02 (2019) 006, arXiv:1810.00007 [hep-ph].
  * [10] M. Hoferichter, B.-L. Hoid, and B. Kubis, “Three-pion contribution to hadronic vacuum polarization,” JHEP 08 (2019) 137, arXiv:1907.01556 [hep-ph].
  * [11] M. Davier, A. Hoecker, B. Malaescu, and Z. Zhang, “A new evaluation of the hadronic vacuum polarisation contributions to the muon anomalous magnetic moment and to $\mathbf{\bm{\alpha}(m_{Z}^{2})}$,” Eur. Phys. J. C 80 no. 3, (2020) 241, arXiv:1908.00921 [hep-ph]. [Erratum: Eur.Phys.J.C 80, 410 (2020)].
  * [12] A. Keshavarzi, D. Nomura, and T. Teubner, “The $g-2$ of charged leptons, $\alpha(M_{Z}^{2})$ and the hyperfine splitting of muonium,” Phys. Rev. D101 (2020) 014029, arXiv:1911.00367 [hep-ph].
  * [13] A. Kurz, T. Liu, P. Marquard, and M. Steinhauser, “Hadronic contribution to the muon anomalous magnetic moment to next-to-next-to-leading order,” Phys. Lett. B734 (2014) 144–147, arXiv:1403.6400 [hep-ph].
  * [14] K. Melnikov and A. Vainshtein, “Hadronic light-by-light scattering contribution to the muon anomalous magnetic moment revisited,” Phys. Rev. D70 (2004) 113006, arXiv:hep-ph/0312226 [hep-ph].
  * [15] P. Masjuan and P. Sánchez-Puertas, “Pseudoscalar-pole contribution to the $(g_{\mu}-2)$: a rational approach,” Phys. Rev. D95 no. 5, (2017) 054026, arXiv:1701.05829 [hep-ph].
  * [16] G. Colangelo, M. Hoferichter, M. Procura, and P. Stoffer, “Dispersion relation for hadronic light-by-light scattering: two-pion contributions,” JHEP 04 (2017) 161, arXiv:1702.07347 [hep-ph].
  * [17] M. Hoferichter, B.-L. Hoid, B. Kubis, S. Leupold, and S. P. Schneider, “Dispersion relation for hadronic light-by-light scattering: pion pole,” JHEP 10 (2018) 141, arXiv:1808.04823 [hep-ph].
  * [18] A. Gérardin, H. B. Meyer, and A. Nyffeler, “Lattice calculation of the pion transition form factor with $N_{f}=2+1$ Wilson quarks,” Phys. Rev. D100 no. 3, (2019) 034520, arXiv:1903.09471 [hep-lat].
  * [19] J. Bijnens, N. Hermansson-Truedsson, and A. Rodríguez-Sánchez, “Short-distance constraints for the HLbL contribution to the muon anomalous magnetic moment,” Phys. Lett. B798 (2019) 134994, arXiv:1908.03331 [hep-ph].
  * [20] G. Colangelo, F. Hagelstein, M. Hoferichter, L. Laub, and P. Stoffer, “Longitudinal short-distance constraints for the hadronic light-by-light contribution to $(g-2)_{\mu}$ with large-$N_{c}$ Regge models,” JHEP 03 (2020) 101, arXiv:1910.13432 [hep-ph].
  * [21] T. Blum, N. Christ, M. Hayakawa, T. Izubuchi, L. Jin, C. Jung, and C. Lehner, “The hadronic light-by-light scattering contribution to the muon anomalous magnetic moment from lattice QCD,” Phys. Rev. Lett. 124 no. 13, (2020) 132002, arXiv:1911.08123 [hep-lat].
  * [22] G. Colangelo, M. Hoferichter, A. Nyffeler, M. Passera, and P. Stoffer, “Remarks on higher-order hadronic corrections to the muon $g-2$,” Phys. Lett. B735 (2014) 90–91, arXiv:1403.7512 [hep-ph].
  * [23] M. Davier, A. Hoecker, B. Malaescu, and Z. Zhang, “Reevaluation of the Hadronic Contributions to the Muon g-2 and to alpha(MZ),” Eur. Phys. J. C 71 (2011) 1515, arXiv:1010.4180 [hep-ph]. [Erratum: Eur.Phys.J.C 72, 1874 (2012)].
  * [24] T. Aoyama et al., “The anomalous magnetic moment of the muon in the Standard Model,” Phys. Rept. 887 (2020) 1–166, arXiv:2006.04822 [hep-ph].
  * [25] G. Bélanger, C. Delaunay, and S. Westhoff, “A Dark Matter Relic From Muon Anomalies,” Phys. Rev. D 92 (2015) 055021, arXiv:1507.06660 [hep-ph].
  * [26] K. Kowalska and E. M. Sessolo, “Expectations for the muon g-2 in simplified models with dark matter,” JHEP 09 (2017) 112, arXiv:1707.00753 [hep-ph].
  * [27] B. Barman, D. Borah, L. Mukherjee, and S. Nandi, “Correlating the anomalous results in $b\to s$ decays with inert Higgs doublet dark matter and muon $(g-2)$,” Phys. Rev. D 100 no. 11, (2019) 115010, arXiv:1808.06639 [hep-ph].
  * [28] L. Calibbi, R. Ziegler, and J. Zupan, “Minimal models for dark matter and the muon g$-$2 anomaly,” JHEP 07 (2018) 046, arXiv:1804.00009 [hep-ph].
  * [29] K. Kowalska and E. M. Sessolo, “Minimal models for g-2 and dark matter confront asymptotic safety,” Phys. Rev. D 103 no. 11, (2021) 115032, arXiv:2012.15200 [hep-ph].
  * [30] S. Jana, P. K. Vishnu, W. Rodejohann, and S. Saad, “Dark matter assisted lepton anomalous magnetic moments and neutrino masses,” Phys. Rev. D 102 no. 7, (2020) 075003, arXiv:2008.02377 [hep-ph].
  * [31] E. J. Chun and T. Mondal, “Explaining $g-2$ anomalies in two Higgs doublet model with vector-like leptons,” JHEP 11 (2020) 077, arXiv:2009.08314 [hep-ph].
  * [32] A. S. De Jesus, S. Kovalenko, F. S. Queiroz, C. Siqueira, and K. Sinha, “Vectorlike leptons and inert scalar triplet: Lepton flavor violation, $g-2$, and collider searches,” Phys. Rev. D 102 no. 3, (2020) 035004, arXiv:2004.01200 [hep-ph].
  * [33] P. Athron, C. Balázs, D. H. Jacob, W. Kotlarski, D. Stöckinger, and H. Stöckinger-Kim, “New physics explanations of $a_{\mu}$ in light of the FNAL muon $g-2$ measurement,” arXiv:2104.03691 [hep-ph].
  * [34] G. Arcadi, L. Calibbi, M. Fedele, and F. Mescia, “Systematic approach to $B$-physics anomalies and $t$-channel dark matter,” arXiv:2103.09835 [hep-ph].
  * [35] X. Qi, A. Yang, W. Liu, and H. Sun, “Scalar dark matter and Muon $g-2$ in a $U(1)_{L_{\mu}-L_{\tau}}$ model,” arXiv:2106.14134 [hep-ph].
  * [36] G. Arcadi, L. Calibbi, M. Fedele, and F. Mescia, “Muon $g-2$ and $B$-anomalies from Dark Matter,” Phys. Rev. Lett. 127 no. 6, (2021) 061802, arXiv:2104.03228 [hep-ph].
  * [37] G. Guedes and J. Santiago, “New leptons with exotic decays: collider limits and dark matter complementarity,” arXiv:2107.03429 [hep-ph].
  * [38] C.-T. Lu, R. Ramos, and Y.-L. S. Tsai, “Shedding light on dark matter with recent muon $(g-2)$ and Higgs exotic decay measurements,” arXiv:2104.04503 [hep-ph].
  * [39] A. E. Cárcamo Hernández, S. Kovalenko, M. Maniatis, and I. Schmidt, “Fermion mass hierarchy and g-2 anomalies in an extended 3HDM Model,” arXiv:2104.07047 [hep-ph].
  * [40] A. E. C. Hernández, S. F. King, and H. Lee, “Fermion mass hierarchies from vectorlike families with an extended 2HDM and a possible explanation for the electron and muon anomalous magnetic moments,” Phys. Rev. D 103 no. 11, (2021) 115024, arXiv:2101.05819 [hep-ph].
  * [41] A. Crivellin and M. Hoferichter, “Consequences of chirally enhanced explanations of $(g-2)_{\mu}$ for $h\to\mu\mu$ and $Z\to\mu\mu$,” JHEP 07 (2021) 135, arXiv:2104.03202 [hep-ph].
  * [42] S. Bhattacharya, P. Poulose, and P. Ghosh, “Multipartite Interacting Scalar Dark Matter in the light of updated LUX data,” JCAP 04 (2017) 043, arXiv:1607.08461 [hep-ph].
  * [43] S. Bhattacharya, P. Ghosh, T. N. Maity, and T. S. Ray, “Mitigating Direct Detection Bounds in Non-minimal Higgs Portal Scalar Dark Matter Models,” JHEP 10 (2017) 088, arXiv:1706.04699 [hep-ph].
  * [44] S. Bhattacharya, P. Ghosh, A. K. Saha, and A. Sil, “Two component dark matter with inert Higgs doublet: neutrino mass, high scale validity and collider searches,” JHEP 03 (2020) 090, arXiv:1905.12583 [hep-ph].
  * [45] S. Profumo, F. Queiroz, and C. Siqueira, “Has AMS-02 Observed Two-Component Dark Matter?,” J. Phys. G 48 no. 1, (2020) 015006, arXiv:1903.07638 [hep-ph].
  * [46] Planck collaboration, N. Aghanim et al., “Planck 2018 results. VI. Cosmological parameters,” (2018) , arXiv:1807.06209 [astro-ph.CO].
  * [47] B. Díaz Sáez, K. Möhling, and D. Stöckinger, “Two Real Scalar WIMP Model in the Assisted Freeze-Out Scenario,” arXiv:2103.17064 [hep-ph].
  * [48] T. Toma, “Internal Bremsstrahlung Signature of Real Scalar Dark Matter and Consistency with Thermal Relic Density,” Phys. Rev. Lett. 111 (2013) 091301, arXiv:1307.6181 [hep-ph].
  * [49] G. Belanger and J.-C. Park, “Assisted freeze-out,” JCAP 03 (2012) 038, arXiv:1112.4491 [hep-ph].
  * [50] G. Belanger, F. Boudjema, and A. Pukhov, “micrOMEGAs : a code for the calculation of Dark Matter properties in generic models of particle interaction,” in Theoretical Advanced Study Institute in Elementary Particle Physics: The Dark Secrets of the Terascale.  2013\.  arXiv:1402.0787 [hep-ph].
  * [51] J. e. a. Aprile, E. Aalbers, “Dark matter search results from a one ton-year exposure of xenon1t,” Physical Review Letters 121 no. 11, (Sep, 2018) . http://dx.doi.org/10.1103/PhysRevLett.121.111302.
  * [52] S. Chang, R. Edezhath, J. Hutchinson, and M. Luty, “Leptophilic Effective WIMPs,” Phys. Rev. D 90 no. 1, (2014) 015011, arXiv:1402.7358 [hep-ph].
  * [53] J. Kawamura, S. Okawa, and Y. Omura, “Current status and muon g-2 explanation of lepton portal dark matter,” JHEP 08 (2020) 042, arXiv:2002.12534 [hep-ph].
  * [54] J. M. Cline, K. Kainulainen, P. Scott, and C. Weniger, “Update on scalar singlet dark matter,” Phys. Rev. D 88 (2013) 055025, arXiv:1306.4710 [hep-ph]. [Erratum: Phys.Rev.D 92, 039906 (2015)].
  * [55] M. Hoferichter, P. Klos, J. Menéndez, and A. Schwenk, “Improved limits for Higgs-portal dark matter from LHC searches,” Phys. Rev. Lett. 119 no. 18, (2017) 181803, arXiv:1708.02245 [hep-ph].
  * [56] T. Bringmann, X. Huang, A. Ibarra, S. Vogl, and C. Weniger, “Fermi LAT Search for Internal Bremsstrahlung Signatures from Dark Matter Annihilation,” JCAP 07 (2012) 054, arXiv:1203.1312 [hep-ph].
  * [57] CMS Collaboration, A. M. Sirunyan et al., “Search for new physics in events with two soft oppositely charged leptons and missing transverse momentum in proton-proton collisions at $\sqrt{s}=$ 13 TeV,” Phys. Lett. B 782 (2018) 440–467, arXiv:1801.01846 [hep-ex].
  * [58] ATLAS Collaboration, G. Aad et al., “A search for the dimuon decay of the Standard Model Higgs boson with the ATLAS detector,” Phys. Lett. B 812 (2021) 135980, arXiv:2007.07830 [hep-ex].
  * [59] CMS Collaboration, A. M. Sirunyan et al., “Evidence for Higgs boson decay to a pair of muons,” JHEP 01 (2021) 148, arXiv:2009.04363 [hep-ex].
  * [60] G. Hiller, C. Hormigos-Feliu, D. F. Litim, and T. Steudtner, “Anomalous magnetic moments from asymptotic safety,” Phys. Rev. D 102 no. 7, (2020) 071901, arXiv:1910.14062 [hep-ph].