11institutetext: Institute of Physics of the Czech Academy of Sciences, Prague
22institutetext: Astronomical Institute of the Czech Academy of Sciences,
Ondřejov 33institutetext: Nicolaus Copernicus Astronomical Center, Polish
Academy of Sciences, Warsaw 44institutetext: Faculty of Mathematics and
Physics, Charles University, Prague

# New developments in aerosol measurements using stellar photometry

Jan Ebr 11<EMAIL_ADDRESS>Jakub Juryšek 11 Michael Prouza 11 Jiří Blažek 11 Petr
Trávníček 11 Dušan Mandát 11 Miroslav Pech for the Pierre Auger Collaboration
and the CTA Consortium Full author list:
http://www.auger.org/archive/authors_2018_09.htmlhttp://www.cta-
observatory.org/11 Petr Janeček for the CTA Consortium33footnotemark: 3 11
Sergey Karpov 11 Ronan Cunniffe 11 Martin Mašek 11 Jiří Eliášek 1144 Martin
Jelínek 22 Ivana Ebrová 1133

###### Abstract

The idea of using stellar photometry for atmospheric monitoring for optical
experiments in high-energy astrophysics is seemingly straightforward, but
reaching high precision of the order of 0.01 in the determination of the
vertical aerosol optical depth (VAOD) has proven difficult. Wide-field
photometry over a large span of altitudes allows a fast determination of VAOD
independently of the absolute calibration of the system, while providing this
calibration as a useful by-product. Using several years of data taken by the
FRAM (F/(Ph)otometric Robotic Atmospheric Monitor) telescope at the Pierre
Auger Observatory in Argentina and about a year of data taken by a similar
instrument deployed at the planned future Southern site of the Cherenkov
Telescope Array in Chile, we have developed methods to improve the precision
of this measurement technique towards and possibly beyond the 0.01 mark.
Detailed laboratory measurements of the response of the whole system to both
the spectrum and intensity of incoming light have proven indispensable in this
analysis as the usual assumption of linearity of the CCD detectors is not
valid anymore for the conditions of the observations.

##  1 Introduction

In the previous AtmoHEAD conference, we have provided a detailed description
VAOD  of the method to measure the vertical aerosol optical depth (VAOD) by
comparing the apparent brightness $m_{\mathrm{inst}}$ (in magnitudes) of a
large amount of stars observed on a series of images from a wide-field
telescope with the expected value $m_{\mathrm{cat}}$ extracted from a catalog.
These images are usually taken to cover a wide range of airmasses $A$ in a
vertical or slanted “scan”. The overall atmospheric extinction is extracted
from a fit of the data with a model that includes both atmospheric and
instrumental effects and can be written, for each star, as:

$\begin{split}m_{\mathrm{inst}}=Mm_{\mathrm{cat}}+Z_{i}+k_{i}A+c_{1}(B-V)(c_{2}(B-V)+1)+\\\
R_{1}r(R_{2}r+1)+k_{C}A(B-V)+k_{A2}A^{2}\end{split}$ (1)

where $B-V$ is the color index of the star, $r$ its distance from the center
of the chip, $M,c_{1},c_{2},R_{1},R_{2},k_{C},k_{A2}$ are global parameters
that depend on the experimental setup, $Z_{i}$ is the calibration constant
(zeropoint) for each scan and $k_{i}$ is the extinction coefficient from which
the VAOD can be deduced after subtracting the effects of the molecular
atmosphere.

Currently, there are two such wide-field telescopes (FRAMs – F/(Ph)otometric
Robotic Atmospheric Monitors) in operation: at the Pierre Auger Observatory
AUGER  near Malargüe, Argentina and at the future planned Southern site of the
Cherenkov Telescope Array (CTA) CTA  near Cerro Paranal, Chile (see PJ  for
details). As the atmospheric conditions on these sites, as well as the
hardware setups, are quite different, comparing the results from the two FRAMs
brings valuable insights for the further development of the method. The
molecular contribution and the $k_{A2}$ parameter must be obtained from the
spectral response of the system which has been, for every combination of
hardware used at any site, measured in the optical laboratory in Olomouc (the
$k_{C}$ parameter can be either calculated or fitted from the data). Note that
all of the measurements presented here are taken in the photometric B filter
at an effective wavelength for VAOD determination between 430 and 440 nm
depending on the type of aerosols and the specific hardware setup.

Figure 1: The dependence of the measured VAOD on the Moon altitude above the
horizon for the Auger (left) and CTA (right) FRAMs. The top row shows the
results before the non-linearity correction, the bottom row shows results with
the non-linearity correction, including the dependence of the bias signal on
outside temperature.

##  2 Moon effect and CCD non-linearity

Since early 2016 the Auger FRAM has started taking data also during the
presence of the Moon on the sky and it quickly became apparent that such data
show a different VAOD distribution from those taken during moonless nights.
This effect has been also confirmed by the CTA FRAM (top row of Fig. 1).
Moreover, a less prominent effect of similar kind is seen when scans are taken
by the Auger FRAM in the direction of the nearby city of Malargüe. This points
to a (spurious) dependence of the results on the sky background – however
attempts to include such an effect in the model have resulted in neither an
improvement of the fit nor the disappearance of the Moon effect.

It turns out that the likely culprit is the non-linearity of small signals of
the CCD cameras we use (Moravian Instruments G4-16000). While this value is
perfectly within the expected performance standards of less than 1% of the
dynamical range (it is less than 0.5 % for any camera we measured), this
translates to up to 30 % of the signal itself for small signals. A typical
observer using such a camera chooses the exposure so that the signal from the
object of interest falls somewhere in the middle of the dynamic range, and
thus experiences a non-linearity of less than 1 % of the signal, whereas in
FRAM operation, where a large number of stars need to be captured in a
reasonable time frame, the signal of many individual stars reaches only a
small fraction of the whole dynamical range, in particular at large airmasses
where extinction is significant. In this range, the behavior of the cameras
was not previously characterized in detail and our laboratory measurements
Sergey  have shown that it is complex and varies from camera to camera.

Figure 2: Comparison between Sun/Moon Photometer and CTA FRAM VAOD
measurements within 15 minutes from each other. Each point represents a FRAM
measurement; several more Photometer measurements usually fall within 15
minutes and their values are averaged. The top row shows the comparison and
histogram of differences for original Photometer data, the bottom row shows
the effect of a fitted correction for Moon phase. The empty gray circles in
the upper left plot show the excluded FRAM outliers that are not considered in
the histogram or the Moon phase fit.

We have thus applied the measured non-linearity corrections (NLC) to all VAOD
data. Because this correction is large and rapidly changes for very small
signal values, care must be taken to properly subtract the bias signal from
each frame, so that only actual light signal is fed to the NLC. This bias
signal depends on the temperature of the readout electronics of the camera,
which is not actively controlled (unlike the CCD chip itself), but is measured
and reported for each image. Using thousands of dark images for each camera,
we have created calibration functions which describe the bias signal
reasonably well; recently we have implemented the readout of masked regions of
the CCD chips of the cameras so that we can measure the bias signal for each
frame directly. Using the NLC without the temperature bias correction renders
the resulting data unsuitable for VAOD analysis altogether. For both Auger and
CTA FRAMs the NLC (with temperature bias correction) removes most of the Moon
effect (Fig. 1) even though the Moon-affected data show a slightly larger
spread. Note that, somewhat counter-intuitively, the NLC mostly changes the
Moonless data, as when Moon is present, the sky background shifts the pixel
levels into higher signal values where the non-linearity is smaller.

The $M$ parameter in the model should have in principle compensated for a
small non-linearity. In reality, the non-linearity at pixel level translates
into a different effect for each individual star depending on the point spread
function (PSF) of the star which is influenced by focus, position in field of
view and any adjustment errors. Interestingly, for CTA FRAM, the apparent non-
linearity determined from the fit ($M-1$) decreases from 2-6 % without the NLC
to <0.5 % with the NLC, while for the Auger FRAM, it is only modestly affected
– this again shows that correcting for this effect at magnitude level is
already too late in the processing chain. It is also important to note that
the non-linearity has been discovered independently (when searching for
variable objects in the data archive) and it cannot be easily inferred from
the VAOD fits; the NLC does not appreciably change the quality of the fit. The
RMS of the spread of individual stars around the fitted model is typically
around 0.1 mag (which is consistent with the published accuracy of the Tycho2
catalog) and thus an effect of the scale of approximately 0.02 mag may not be
apparent. Thus a naive determination of the VAOD from the fit may lead to a
wrong conclusion even though no obvious discrepancy (besides the Moon effect)
is observed.

##  3 Precision of the measurements

The statistical error in the determination of the VAOD from an individual scan
is between 0.003 and 0.008 for the Auger FRAM and between 0.002 and 0.004 for
the CTA FRAM which can, due to its larger field of view, cover a larger span
in airmass in a given time window. The precision is thus largely determined by
systematics. Several sources are easily quantified: the subtraction of the
molecular atmosphere (0.003, shall be improved using the MODTRAN software and
GDAS atmospheric profiles GDAS ), the uncertainty in determining the system
parameters from the fit (0.003), the freedom of choice of the cutoffs in
brightness and airmass (0.005). Further sources are yet to be quantified: the
precision of the measurements of the system spectral response, the quality of
the bias correction, possible biases in the stellar catalog, the influence of
the automatic rejection of outlying stars, residual cloud contamination, etc.
The effect causing the systematic outliers shown in PJ  also needs to be
identified and either removed or included in the systematics (unless it is
physical, which is unlikely, see next section). Any possible significant
improvement in the systematics will inevitably mean that the statistical error
becomes dominant – this is however to a large extent likely a feature of the
Tycho2 catalog and the attempts to use different catalogs have been met with
mixed results so far due to coverage inhomogenity (APASS) or very different
spectral bands (GAIA) and thus this might be the ultimate limit of the method
until a better all-sky photometric catalog is available.

##  4 Sun/Moon photometer comparison

Comparison of the FRAM VAOD data with other instruments, such as LIDARs, is
notoriously difficult because the FRAM measures the integral extinction all
the way through the atmosphere. However the Sun/Moon Photometer Cimel CE318-T,
installed next to the CTA FRAM in Chile as a part of the same site-
chracterisation campaign, does measure the same quantity using the Sun and the
Moon Jakub  and it is thus ideal for this comparison. Unfortunately, due to
communication problems on the remote location, the data from this instrument
could be retrieved only after a recent on-site visit and thus the data
processing presented here is extremely preliminary.

Fig. 2 shows the comparison between FRAM and Photometer measurements. The
Photometer measures VAOD in many narrow-band filters from which we have chosen
the 440 nm one which is very close (within 10 nm) to the effective mean
wavelength of all FRAM setups. The outliers in the FRAM data shown in PJ  were
removed for this purpose. The mean value for
$VAOD_{\mathrm{FRAM}}-VAOD_{\mathrm{Photometer}}$ is 0.009 with 68 % of the
data contained within $\pm 0.02$, which is highly encouraging since the
possible uncertainty of the Moon measurements of the Photometer is quoted as
0.04 (in contrast to <0.01 for Sun measurements). Nevertheless, an obvious
dependence on the Moon phase is visible in the data (across several lunar
cycles), pointing to a small deficiency in the Moon illumination model. We
have fitted two separate fourth-order polynomials to the difference between
the FRAM and Photometer data for waxing and gibbous Moon (as there is a large
step in this difference around full Moon). Applying this correction reduces
scatter of the data and 68 % of the differences is now contained within $\pm
0.016$ (with the mean value being 0 by construction, see bottom row of Fig.
2).

The Sun/Moon Photometer was also taken for three months in 2017 to Argentina,
where it performed a dedicated measurement campaign next to the Auger FRAM.
Unfortunately, the weather in this period was less than optimal and the amount
of useful data taken is small. The data also show an apparent dependence on
lunar phase, but because for many phases only one data point was obtained, it
is not possible to do a similar fit. We cannot use the correction obtained
from the CTA FRAM as the Photometer is calibrated for each site individually
using the respective data.

## Acknowledgments

The authors are grateful for the support by the grants of the Ministry of
Education, Youth and Sports of the Czech Republic, projects LTT17006,
LM2015046, LTT18004, LM2015038, CZ.02.1.01/0.0/0.0/16_013/0001402, CZ.
02.1.01/0.0/0.0/16_013/0001403 and CZ.02.1.01/0.0/0.0 /15_003/0000437. This
research was supported in part by the Polish National Science Centre under
grant 2017/26/D/ST9/00449. The work is partially performed according to the
Russian Government Program of Competitive Growth of Kazan Federal University.

## References

  * (1) Ebr, J. et al. (The Pierre Auger Collaboration), EPJ Web Conf. 144 01011 (2017)
  * (2) Janecek, P. et al. (The Pierre Auger Collaboration and CTA Consortium), these proceedings
  * (3) The Pierre Auger Collaboration, Nucl. Instrum. Meth. A 798, 172 (2015)
  * (4) CTA Consortium, Astroparticle Physics 43, pp. 3–18 (2013)
  * (5) Karpov, S. V. et al., Astronomische Nachrichten, 2018, in press
  * (6) Ebr, J. et al., Proceedings of the 3rd Astrorob Workshop, Rev. Mex. AA (SC) 45, 114 (2014)
  * (7) The Pierre Auger Collaboration, Journal of Instrumentation 7, P09001 (2012)
  * (8) Ebr, J. et al. (The Pierre Auger Collaboration), Rev. Mex. AA (SC) 45, 53 (2014)
  * (9) Munar, P., Gaug, M. (CTA Consortium), these proceedings
  * (10) Jurysek, J. et al (CTA Consortium) PoS ICRC2017 (2018) 853


# Andreev bound states in Superconductor-Quantum dot Josephson junction at
infinite-U limit

Tanuj Chamoli<EMAIL_ADDRESS>Ajay<EMAIL_ADDRESS>Department
of Physics, Indian Institute of Technology Roorkee : 247667, Uttarakhand,
India

###### Abstract

Andreev bound states (ABSs) are studied in quantum dot coupled to conventional
BCS superconducting leads on the basis of effective slave boson Hamiltonian in
an infinite-U (Coulomb interaction) limit followed by Green’s function
technique. From the relevant Green’s function, density of states (DOS) is
analyzed at different superconducting phase differences. On the basis of
numerical computation, it is pointed out from DOS plot that ABSs (sub-gap
states) arise due to the phase difference between left and right
superconducting leads. Due to the dependence on superconducting phase
difference, ABSs are current carrying states. We have also analyzed the energy
of ABSs as a function of phase difference between left and right
superconducting lead for electron hole symmetric as well as non symmetric
cases at various dot-lead coupling strengths. It is also pointed out that a
finite internal gap arises between upper and lower Andreev bound states in the
absence of electron-hole or left-right symmetry. These results are viewed in
the light of existing theoretical analysis.

###### keywords:

, Quantum dot, BCS Superconductor, Andreev bound states, Infinite-U slave
Boson mean field approach

###### MSC:

[2020] 81Q37, 81T99 , 81V65 , 82D55

††journal: Superlattices and Microstructures

##  1 Introduction

Josephson transport across quantum dot junction has been widely used in
nanoelectronic devices such as in superconducting quantum interference devices
(SQUIDS) [1], Cooper pair splitters [2] and superconducting quantum bits [3]
which are the basis for superconducting quantum computers. Josephson junction
is obtained by sandwiching thin insulating layer between two superconducting
leads and Josephson supercurrent flows across the junction due to the phase
difference between Cooper pairs present on the left and right superconducting
lead [4]. Insulating layer can be replaced by some nanoscopic structures such
as quantum dots (QDs) [5]. QD itself, is not a superconducting but Cooper
pairs tunneling can take place through it because coherence length of Cooper
pairs is significantly larger as compared to the size of the QD. Due to the
proximity effect [6], induced pairing takes place on the dot around the Fermi
level. In this way, a reduced gap is arised in the dot. Consequently, subgap
states are observed in S-QD-S system which actually arises due to the
superconducting phase difference between two leads [7, 8]. These states are
called Andreev bound states and are responsible for supercurrent flowing
through the junction. Andreev bound states are consequences of the phenomena
called multiple Andreev reflection [9].

For odd number of electrons occupation , QD acts as a magnetic impurity, and
spin of quantum dot gets coupled to spin of one electron of the Cooper pair
from the leads with a characteristic binding energy $k_{B}T_{k}$, where
$T_{K}$ is the Kondo temperature. Kondo temperature for normal metallic lead
varies as $T_{k}\sim e^{-1/{\rho J}}$, where $J$ is the exchange interaction
between spin of QD and spin of one of the electron of metallic lead [10]. It
leads to the formation of Kondo singlet in N-QD-N junction and it can be
visualized as a peak in density of states. Although for superconducting leads,
there is an attaractive interaction between two electrons to form a bound
Cooper pair. Therefore, interplay arises between multiple Andreev reflection
and Kondo effect in S-QD-S junction [11, 12, 13]. It is observed
experimentally that for $T_{k}>\Delta$, where $\Delta$ is binding energy of
Cooper pair, there is always likelihood of formation of ABSs in the dot. It
give rise to sharp increase in DOS within superconducting gap at $+\omega_{s}$
and $-\omega_{s}$, where $\omega_{s}$ corresponds to energy of Andreev bound
states. Therefore, the ratio $\Delta/T_{k}$ is an important parameter for the
transmission of Cooper pair from left to right lead [14]. For even number of
electrons occupation on dot, positive supercurrent while for odd number of
electrons occupation negative supercurrent flows through the junction. It
indicates that the supercurrent reverses the sign by adding an electron on the
QD [15].

Figure 1: Schematic of S-QD-S junction, showing coupling of QD with
conventional superconducting leads with respect to Fermi level (shown as red
dashed line)

As of now, several experimental attempts have been made to fabricate and
characterize nanowire devices in which nanostructure is coupled to
superconducting leads [16, 17, 18, 19]. These experiments indicate the
existence of subgap states in S-QD-S system at various parameter
regime.Further it is analyzed that these subgap states are effected by
superconducting phase difference and dot-lead coupling strength. Outcomes of
these experiments are analyzed by various analytical and computational
techniques, which are based on different approximations. To mention few of
these theoretical techniques include renormalization group approach [20, 21,
22], quantum monte carlo simulation [23, 24], non-crossing approximation [25,
26, 27], perturbation theory [28] and also non-equlibrium Green’s function
[29, 30]. Although infinite-U slave boson mean field (SBMF) method to
understand the influence of infinite Coulomb interaction on the QD state
connected to BCS superconducting leads is addressed by few of the researchers
[31, 32], but as of now , role of strong onsite interaction on the formation
of subgap states and on nature of Josephson transport as well is not clearly
understood from theoretical point of view. Also extension of S-QD-S junction
to multiple coupled QDs is yet to be analyzed. Therefore, in present attempt
we have planned to address the issue related to strong on site Coulomb
interaction (infinite-U) on the ABSs and their phase dependence, where we have
used Bogolibouv transformation followed by infinite-U slave Boson mean field
approach to understand the nature of S-QD-S junction at infinite Coulomb
interaction. In the preceding section, we have presented the theoretical model
for S-QD-S system.

##  2 Theoretical formulation

As depicted in Fig. 1 the model Hamiltonian for S-QD-S junction where QD has
single level is modelled by single impurity Anderson Hamiltonian (SIAM). BCS
term is introduced to take care of superconducting leads.

$H=H_{leads}+H_{QD}+H_{dot-lead}$ (1)

where contribution of leads is given by

$H_{leads}=\sum_{\alpha k\sigma}{\epsilon_{\alpha
k\sigma}{c^{\dagger}}_{\alpha k\sigma}c_{\alpha k\sigma}}-\sum_{\alpha
k\sigma}({\Delta_{k\alpha}{c^{\dagger}}_{\alpha
k\sigma}{c^{\dagger}}_{\alpha-k-\sigma}+h.c.)}$ (2)

where $\epsilon_{k\alpha\sigma}$ denotes the energy of an electron with spin
$\sigma$ ($\sigma=\uparrow,\downarrow$ denotes spin of spin up and spin down
electrons respectively) in the $k^{\text{th}}$ energy state of lead $\alpha$
($\alpha=L,R$ corresponds to left and right lead).
$\Delta_{k\alpha}=\left|\Delta_{k\alpha}\right|e^{i\phi_{k\alpha}}$ is complex
superconducting order parameter in which $\phi_{k\alpha}$ is phase
correponding to lead $\alpha$. For BCS, S-wave superconductor, superconducting
order parameter does not have k dependency. Quantum Dot contribution is given
as follows.

$H_{QD}=\sum_{\sigma}\epsilon_{d}{d^{\dagger}}_{\sigma}d_{\sigma}$ (3)

$\epsilon_{d}$ is the energy level of the quantum dot. Coupling between dot
level and superconducting lead is of the form

$H_{dot-lead}=\sum_{\alpha k\sigma}(h_{k\alpha}{c^{\dagger}}_{\alpha
k\sigma}d_{\sigma}+h_{k\alpha}^{*}{d{{}^{\dagger}}}_{\sigma}c_{\alpha
k\sigma})$ (4)

where $h_{k\alpha}$ is coupling strength between dot and lead. To diagonalize
the Hamiltonian Bogoliubov transformation is used, defined as follows

$c_{k\alpha\uparrow}=u_{k\alpha}^{*}\beta_{{k\alpha}\uparrow}+v_{k\alpha}{\beta^{\dagger}}_{-{k\alpha}\downarrow}$
(5)

$c_{-k\alpha\downarrow}=u_{k\alpha}^{*}\beta_{-{k\alpha}\downarrow}-v_{k\alpha}{\beta^{\dagger}}_{{k\alpha}\uparrow}$
(6)

${c^{\dagger}}_{k\alpha\uparrow}=u_{k\alpha}{\beta^{\dagger}}_{{k\alpha}\uparrow}+v_{k\alpha}^{*}{\beta}_{-{k\alpha}\downarrow}$
(7)

${c^{\dagger}}_{-k\alpha\downarrow}=u_{k\alpha}{\beta^{\dagger}}_{-{k\alpha}\downarrow}-v_{k\alpha}^{*}{\beta}_{{k\alpha}\uparrow}$
(8)

where

$\begin{split}{|u_{k\alpha}|}^{2}=\frac{1}{2}\bigg{(}1+\frac{\epsilon_{k\alpha}}{\sqrt{{\epsilon_{k\alpha}}^{2}+{|\Delta_{k\alpha}|}^{2}}}\bigg{)}\\\
{|v_{k\alpha}|}^{2}=\frac{1}{2}\bigg{(}1-\frac{\epsilon_{k\alpha}}{\sqrt{{\epsilon_{k\alpha}}^{2}+{|\Delta_{k\alpha}|}^{2}}}\bigg{)}\\\
u_{k\alpha}^{*}v_{k\alpha}=\frac{\Delta_{k\alpha}}{2E_{k\alpha}}\\\
v_{k\alpha}^{*}u_{k\alpha}=\frac{\Delta_{k\alpha}^{*}}{2E_{k\alpha}}\end{split}$
(9)

After Bogoliubov transformation, Hamiltonian takes the form as follows.

$H_{Bog}=\sum_{\alpha k}{E_{\alpha k}({\beta^{\dagger}}_{\alpha
k\uparrow}\beta_{\alpha
k\uparrow}}+{\beta^{\dagger}}_{\alpha-k\downarrow}\beta_{\alpha-k\downarrow})+\sum_{\sigma}\epsilon_{d}{d^{\dagger}}_{\sigma}d_{\sigma}\\\
+\sum_{\alpha
k}[h_{k\alpha}(u_{k\alpha}{\beta^{\dagger}}_{k\alpha\uparrow}+v_{k\alpha}^{*}{\beta}_{-k\alpha\downarrow})d_{\uparrow}+(u_{k\alpha}{\beta^{\dagger}}_{k\alpha\downarrow}-v_{k\alpha}^{*}{\beta}_{-k\alpha\uparrow})d_{\downarrow}]\\\
+h_{k\alpha}^{*}[d^{\dagger}_{\uparrow}(u_{k\alpha}^{*}\beta_{k\alpha\uparrow}+v_{k\alpha}{\beta^{\dagger}}_{-k\alpha\downarrow})+d^{\dagger}_{\downarrow}(u_{k\alpha}^{*}\beta_{k\alpha\downarrow}-v_{k\alpha}{\beta^{\dagger}}_{-k\alpha\uparrow})]$
(10)

where

$E_{k\alpha}={\sqrt{{\epsilon_{k\alpha}}^{2}+{|\Delta_{k\alpha}|}^{2}}}$ (11)

Now to study the system at infiite-U employing slave boson mean field
approach, auxiliary fermionic operators $f_{\sigma}$ and
${f^{\dagger}}_{\sigma}$ are introduced in terms of auxiliary boson field
operators $b$ and ${b^{\dagger}}$.

${d_{\sigma}}={b^{\dagger}}f_{\sigma}$ (12)

${d^{\dagger}}_{\sigma}={f^{\dagger}}_{\sigma}b$ (13)

Due to infinite-U, single electron occupancy will be there on QD. Condition
for single occupancy is given as

${b^{\dagger}}b+\sum_{\sigma}{f^{\dagger}}_{\sigma}f_{\sigma}=1$ (14)

This condition allows to apply the slave-boson technique. Hamiltonian now
takes the form as

$H^{{}^{\prime}}=\sum_{\alpha k}{E_{\alpha k}({\beta^{\dagger}}_{\alpha
k\uparrow}\beta_{\alpha
k\uparrow}}+{\beta^{\dagger}}_{\alpha-k\downarrow}\beta_{\alpha-k\downarrow})+\sum_{\sigma}\epsilon_{d}{f^{\dagger}}_{\sigma}bb^{\dagger}f_{\sigma}\\\
+\sum_{\alpha
k}h_{k\alpha}[(u_{k}{\beta^{\dagger}}_{k\uparrow}+{\beta}_{-k\downarrow})f_{\uparrow}b^{\dagger}+{f^{\dagger}}_{\uparrow}(u_{k}\beta_{k\uparrow}+v_{k}{\beta^{\dagger}}_{-k\downarrow})b]\\\
+h_{k\alpha}^{*}[(u_{k}{\beta^{\dagger}}_{k\downarrow}-v_{k}{\beta}_{-k\uparrow})f_{\downarrow}b^{\dagger}+{f^{\dagger}}_{\downarrow}(u_{k}\beta_{k\downarrow}-v_{k}{\beta^{\dagger}}_{-k\uparrow})b]+\lambda\sum_{\sigma}({f^{\dagger}}_{\sigma}f_{\sigma}+b^{\dagger}b-1)$
(15)

where $\lambda$ is Lagrange multiplier. In mean field approximation we replace
Boson operators by their expectation values, $\langle
b^{\dagger}(t)\rangle=\langle b(t)\rangle=b$. Effective Slave Boson mean field
Hamiltonian finally comes out to be as

$H_{SBMF}=\sum_{\alpha k}{E_{\alpha k}({\beta^{\dagger}}_{\alpha
k\uparrow}\beta_{\alpha
k\uparrow}}+{\beta^{\dagger}}_{\alpha-k\downarrow}\beta_{\alpha-k\downarrow})+\sum_{\sigma}\tilde{\epsilon}_{d}{f^{\dagger}}_{\sigma}f_{\sigma}\\\
+\sum_{\alpha
k}\tilde{h_{k\alpha}}[(u_{k}{\beta^{\dagger}}_{k\uparrow}+{\beta}_{-k\downarrow})f_{\uparrow}+{f^{\dagger}}_{\uparrow}(u_{k}\beta_{k\uparrow}+v_{k}{\beta^{\dagger}}_{-k\downarrow})]\\\
+\tilde{h_{k\alpha}^{*}}[(u_{k}{\beta^{\dagger}}_{k\downarrow}-v_{k}{\beta}_{-k\uparrow})f_{\downarrow}+{f^{\dagger}}_{\downarrow}(u_{k}\beta_{k\downarrow}-v_{k}{\beta^{\dagger}}_{-k\uparrow})]+\lambda(b^{2}-1)$
(16)

where $\tilde{\epsilon_{d}}(=\epsilon_{d}+\lambda)$ is renormalized QD level
and $\tilde{h_{K\alpha}}(=bh_{k\alpha})$ is renormalized dot lead coupling
strength. Now, employing Green’s function equation of motion, set of coupled
equations is obtained.

$(\omega-\tilde{\epsilon}_{d})\langle\langle
f_{\uparrow},{f^{\dagger}}_{\uparrow}\rangle\rangle=1+\sum_{k\alpha}\tilde{h_{k\alpha}}^{*}{u_{k\alpha}}^{*}\langle\langle\beta_{k\alpha\uparrow},{f^{\dagger}}_{\uparrow}\rangle\rangle+\sum_{k\alpha}\tilde{h_{k\alpha}}^{*}v_{k\alpha}\langle\langle\beta^{\dagger}_{-k\alpha\downarrow},{f^{\dagger}}_{\uparrow}\rangle\rangle$
(17)

$(\omega-
E_{k\alpha})\langle\langle\beta_{k\alpha\uparrow},{f^{\dagger}}_{\uparrow}\rangle\rangle=\tilde{h_{k\alpha}}u_{k\alpha}\langle\langle
f_{\uparrow},{f^{\dagger}}_{\uparrow}\rangle\rangle+\tilde{h_{k\alpha}}^{*}v_{k\alpha}\langle\langle
f^{\dagger}_{\downarrow},{f^{\dagger}}_{\uparrow}\rangle\rangle$ (18)

$(\omega+E_{k\alpha})\langle\langle\beta^{\dagger}_{-k\alpha\downarrow},{f^{\dagger}}_{\uparrow}\rangle\rangle=\tilde{h_{k\alpha}}{v_{k\alpha}}^{*}\langle\langle
f_{\uparrow},{f^{\dagger}}_{\uparrow}\rangle\rangle-\tilde{h_{k\alpha}}^{*}{u_{k\alpha}}^{*}\langle\langle
f^{\dagger}_{\downarrow},{f^{\dagger}}_{\uparrow}\rangle\rangle$ (19)

$(\omega+\tilde{\epsilon_{d}})\langle\langle
f^{\dagger}_{\downarrow},{f^{\dagger}}_{\uparrow}\rangle\rangle=\sum_{k\alpha}\tilde{h_{k\alpha}}{v_{k\alpha}}^{*}\langle\langle\beta_{-k\alpha\uparrow},{f^{\dagger}}_{\uparrow}\rangle\rangle-\sum_{k\alpha}\tilde{h_{k\alpha}}u_{k\alpha}\langle\langle\beta^{\dagger}_{k\alpha\downarrow},{f^{\dagger}}_{\uparrow}\rangle\rangle$
(20)

$(\omega-
E_{k\alpha})\langle\langle\beta_{-k\alpha\uparrow},{f^{\dagger}}_{\uparrow}\rangle\rangle=\tilde{h_{k\alpha}}u_{k\alpha}\langle\langle
f_{\uparrow},{f^{\dagger}}_{\uparrow}\rangle\rangle+\tilde{h_{k\alpha}}^{*}v_{k\alpha}\langle\langle
f^{\dagger}_{\downarrow},{f^{\dagger}}_{\uparrow}\rangle\rangle$ (21)

$(\omega+E_{k\alpha})\langle\langle\beta^{\dagger}_{k\alpha\downarrow},{f^{\dagger}}_{\uparrow}\rangle\rangle=\tilde{h_{k\alpha}}{v_{k\alpha}}^{*}\langle\langle
f_{\uparrow},{f^{\dagger}}_{\uparrow}\rangle\rangle-\tilde{h_{k\alpha}}^{*}{u_{k\alpha}}^{*}\langle\langle
f^{\dagger}_{\downarrow},{f^{\dagger}}_{\uparrow}\rangle\rangle$ (22)

Solving above coupled equations from equation (17) to (22) , relevant Grren’s
function $\langle\langle f_{\uparrow},{f^{\dagger}}_{\uparrow}\rangle\rangle$
is obtained as follows.

$\langle\langle
f_{\uparrow},{f^{\dagger}}_{\uparrow}\rangle\rangle=\frac{1}{\omega-\tilde{\epsilon}_{d}-\tilde{\Gamma}_{2L}-\tilde{\Gamma}_{2R}-\frac{(\tilde{\Gamma}_{3L}+\tilde{\Gamma}_{3R})(\tilde{\Gamma}_{4L}+\tilde{\Gamma}_{4R})}{\omega+\tilde{\epsilon}_{d}-\tilde{\Gamma}_{1L}-\tilde{\Gamma}_{1R}}}$
(23)

where

$\tilde{\Gamma}_{1\alpha}=\sum_{k}|\tilde{h_{k\alpha}}|^{2}\left[\frac{|v_{k\alpha}|^{2}}{\omega-
E_{k\alpha}}+\frac{|u_{k\alpha}|^{2}}{\omega+E_{k\alpha}}\right]$ (24)

$\tilde{\Gamma}_{2\alpha}=\sum_{k}|\tilde{h_{k\alpha}}|^{2}\left[\frac{|u_{k\alpha}|^{2}}{\omega-
E_{k\alpha}}+\frac{|v_{k\alpha}|^{2}}{\omega+E_{k\alpha}}\right]$ (25)

$\tilde{\Gamma}_{3\alpha}=\sum_{k}{\tilde{h_{k\alpha}}}^{2}v_{k\alpha}^{*}u_{k\alpha}\left[\frac{1}{\omega-
E_{k\alpha}}-\frac{1}{\omega+E_{k\alpha}}\right]$ (26)

$\tilde{\Gamma}_{4\alpha}=\sum_{k}{\tilde{h_{k\alpha}^{*}}}^{2}u_{k\alpha}^{*}v_{k\alpha}\left[\frac{1}{\omega-
E_{k\alpha}}-\frac{1}{\omega+E_{k\alpha}}\right]$ (27)

To obtain the parameters b and $\lambda$, we followed the minimization of free
energy function, which can be described in terms of density of states
[$\tilde{\rho}_{f}(\omega)$] as follows.

$F=-\frac{1}{\beta}lnZ=-\frac{N}{\beta}\int_{-\infty}^{\infty}ln(1+e^{-\beta(\omega-\mu)})\tilde{\rho}_{f}(\omega)d\omega$
(28)

For a wide flat conduction band of width 2D, change in free energy due to
impurity is given as follows [33].

$F_{MF}=-\frac{N}{\beta}\int_{-D}^{D}f(\omega)tan^{-1}(\frac{\tilde{\Lambda}}{\tilde{\epsilon_{d}}-\omega})d\omega+\lambda(b^{2}-1)$
(29)

where $\tilde{\Lambda}$ is renormalized hybridization given as follows

$\tilde{\Lambda}=\tilde{\Gamma}_{2L}+\tilde{\Gamma}_{2R}+\frac{(\tilde{\Gamma}_{3L}+\tilde{\Gamma}_{3R})(\tilde{\Gamma}_{4L}+\tilde{\Gamma}_{4R})}{\omega+\tilde{\epsilon}_{d}-\tilde{\Gamma}_{1L}-\tilde{\Gamma}_{1R}}$
(30)

For flat conduction band, dot-lead coupling becomes k-independent. It is clear
that for b=1 and $\lambda=0$ , $\tilde{\Lambda}$ reduces to that of non
interacting hybridization. Minimization of free energy with respect to
parameters b and $\lambda$ , give rise to two coupled integral equations
described as follows.

$\frac{\partial
F_{MF}}{\partial\lambda}=\frac{N}{\pi}\int_{-D}^{D}\frac{f(\omega)d\omega}{{(\tilde{\epsilon_{d}}-\omega})^{2}+{\tilde{\Lambda}}^{2}}\left[\tilde{\Lambda}+\frac{(\tilde{\Gamma}_{3L}+\tilde{\Gamma}_{3R})(\tilde{\Gamma}_{4L}+\tilde{\Gamma}_{4R})(\tilde{\epsilon_{d}}-\omega)}{(\omega+\tilde{\epsilon}_{d}-\Gamma_{1L}-\tilde{\Gamma}_{1R})^{2}}\right]\\\
+(b^{2}-1)=0$ (31) $\frac{\partial F_{MF}}{\partial
b}=-\frac{N}{\pi}\int_{-D}^{D}\frac{f(\omega)d\omega}{{(\tilde{\epsilon_{d}}-\omega})^{2}}[(\Gamma_{2L}+\Gamma_{2R})+\frac{(\Gamma_{3L}+\Gamma_{3R})(\Gamma_{4L}+\Gamma_{4R})}{(\omega+\tilde{\epsilon_{d}}-\tilde{\Gamma_{1L}}-\tilde{\Gamma_{1R}})^{2}}\\\
(\lambda+\Gamma_{11L}+\Gamma_{1R}+2b(\omega+\tilde{\epsilon_{d}}-\tilde{\Gamma_{1L}}-\tilde{\Gamma_{1R}}))]+(b^{2}-1)=0$
(32)

In equations (32) and (LABEL:fminlam), N is level degenracy ( N=2, for single
level case) and $f(\omega)$ is distribution function. The Bogoliubonos follow
the Fermi-Dirac distribution function for dispersion $E_{k}$, given as
follows.

$f(\omega)=\langle\beta^{\dagger}_{k\uparrow}\beta_{k\uparrow}\rangle=\langle\beta^{\dagger}_{-k\downarrow}\beta_{-k\downarrow}\rangle=\frac{1}{e^{\beta
E_{k}}+1}$ (33)

We have used zero temperature for our calculation. To analyse the DOS and
energy of ABSs, we need to solve equations (32) and (LABEL:fminlam)
numerically for various parameters of SBMF Hamiltonian.

##  3 Results and Discussions

On the basis of SIAM, we have used Bogoliubov transformation followed by
infinite-U slave Boson mean field approach to get the relevant Green’s
function. From Green’s function, density of states of quantum dot is
calculated as follows.

$\rho_{f}(\omega)=-\frac{1}{\pi}imag(G_{f}(\omega^{+}))$ (34)

where, Green’s function $G_{f}(\omega^{+})$ is given by equation (23).

(a) Normalized density of states in QD

(b) close up view of normalized DOS in QD

Figure 2: Normalized DOS (Density of states) as a function of $\omega/\Delta$
for electron-hole symmetric ($\epsilon_{d}$ = 0) as well as left-right
symmetric ($h_{L}=h_{R}$) case at different values of superconducting phase
difference ($\phi=0,\pi/4,\pi/2,3\pi/4,\pi$).

In Fig. (2), normalized density of states (DOS) is analyzed for identical left
and right superconducting lead
($\left|\Delta_{L}\right|=\left|\Delta_{R}\right|=1meV$) as a function of
$\omega/\Delta$ for various values of superconducting phase difference
($\phi=\phi_{L}-\phi_{R}=0,\pi/4,\pi/2,3\pi/4,\pi$) in particle-hole symmetric
($\epsilon_{d}=0$) and left-right symmetric ($h_{L}=h_{R}$) case. It is
observed that in the absence of phase difference between left and right lead
($\phi=0$), superconducting gap can be seen in DOS plot with no subgap states
around the Fermi level. As the value of phase difference is increased to
$\phi=\pi/4$, two sub gap states (symmetrically above and below the Fermi
level) appear in the gap, which are called Andreev bound states. Sub gap state
below the Fermi level is called lower Andreev bound state and one above the
Fermi level is called upper Andreev bound state. In the other words, phase
difference introduces paired state in the quantum dot with reduced gap.
Further increasing the value of $\phi$ ($\phi=\pi/2,3\pi/4$) upper and lower
ABSs come closer and closer. For $\phi=\pi$ these two states merge into each
other and only a single subgap state is obtained. Fig. (2) shows the close up
look of Fig. (2). Nature of DOS can be observed clearly in Fig. (2) at
superconducting gap edge. Here, it can be pointed out that the DOS plot for
infinite-U case is similar to non interacting case except the amplitude of DOS
is decreased. It denotes that transmission of Cooper pair electrons from
superconducting lead to QD is decreased. Hence, ABSs are not effected by
infinite Coulomb interaction in the dot and still present within the
superconducting gap. Dependency of ABSs on phase difference indicates that the
transmission of Cooper pairs can be controlled by changing superconducting
phase difference between left and right lead. These results are valid in the
region $T_{K}>\Delta$, where $T_{K}$ is the Kondo temperature corresponding to
superconducting lead.

Further, the energy of ABSs is calculated from poles of the Green’s function
(23), which is given by the solution of non linear equation.

$(\omega_{s}-\tilde{\epsilon}_{d}-\tilde{\Gamma}_{2L}-\tilde{\Gamma}_{2R})(\omega_{s}+\tilde{\epsilon}_{d}-\tilde{\Gamma}_{1L}-\tilde{\Gamma}_{1R})=(\tilde{\Gamma}_{3L}+\tilde{\Gamma}_{3R})(\tilde{\Gamma}_{4L}+\tilde{\Gamma}_{4R})$
(35)

(a) electron-hole symmetric case

(b) electron-hole non symmetric case

Figure 3: ABS as a function of phase difference at different dot-lead coupling
ratios ($h_{L}/h_{R}=1,1.5,2,2.5$).

Equation (35) is solved numerically. In Fig (3), energy of Andreev bound
states is analyzed with variation in superconducting phase difference in
electron-hole symmetric case for various left lead-dot coupling ($h_{L}$) to
right lead-dot couling ($h_{R}$) ratios. For $h_{L}=h_{R}$, i.e. left-right
symmetry, junction acts as a perfectly transmitting channel. Further
increasing the ratio $h_{L}/h_{R}$ (1.5,2,2.5) a finite internal gap arises
between upper and lower ABSs at $\phi=\pi$. This internal gap increases with
increase in the ratio $h_{L}/h_{R}$. Fig (3) shows the energy of ABSs with
variation in superconducting phase difference in the absence of electron-hole
symmetry for various left lead-dot coupling($h_{L}$) to right lead-dot
coupling ($h_{R}$) ratios. In this case, there is an small internal gap is
present, even at $h_{L}=h_{R}$. Further increasing the ratio $h_{L}/h_{R}$
(1.5,2,2.5) width of this finite internal gap is increased. So,it can be
pointed out that non-perfect transmitting channel is formed at the junction in
the absence of electron-hole symmetry. Due to dependency of ABSs on
superconducting phase difference, these are current carrying states. Hence,
ABS forms a channel through which Josephson current flows across the S-QD-S
junction even at infinite on dot Coulomb interaction limit. These results are
in accordance with the works pointed out in other theoretical works, where sub
gap states in S-QD-S system is studied [32].

##  4 Conclusions and outlook

Employing the infinite-U slave Boson mean field approach for Single impurity
Anderson model of superconductor representing S-QD-S junction, sub gap states
(Andreev bound states) are observed within the superconducting gap. These sub
gap states are found to depend on the superconducting phase difference between
two leads as well as left lead-dot coupling to right lead-dot coupling ratio.
Due to dependence on phase difference, these states are current carrying
states. This formalism can be extended to multilevel quantum dot and also to
multi coupled QDs connected to superconducting leads. More, works is required
for these types of system because of their potential application in nano
electronic devices and hence we are attempting the works in this direction.

## Aknowledgement

One of the authors, Tanuj Chamoli would like to thank Ministry of Human
Resource Development (MHRD), Government of India for financial support,
received as a research fellowship.

## References

  * [1] C.-P. Yang, S.-I. Chu, S. Han, Possible realization of entanglement, logical gates, and quantum-information transfer with superconducting-quantum-interference-device qubits in cavity qed, Phys. Rev. A 67 (2003) 042311.  doi:10.1103/PhysRevA.67.042311.
URL https://link.aps.org/doi/10.1103/PhysRevA.67.042311

  * [2] L. G. Herrmann, F. Portier, P. Roche, A. L. Yeyati, T. Kontos, C. Strunk, Carbon nanotubes as cooper-pair beam splitters, Phys. Rev. Lett. 104 (2010) 026801.  doi:10.1103/PhysRevLett.104.026801.
URL https://link.aps.org/doi/10.1103/PhysRevLett.104.026801

  * [3] W. F. Clarke, J., Superconducting quantum bits, Nature 453 (1) (2008) 1031.
  * [4] B. Josephson, Coupled superconductors, Reviews of Modern Physics 36 (1) (1964) 216\.
  * [5] M. A. Reed, Quantum dots, Scientific American 268 (1) (1993) 118–123.
  * [6] Y. Avishai, A. Golub, A. D. Zaikin, Tunneling through an anderson impurity between superconductors, Physical Review B 63 (13) (2001) 134515.
  * [7] Y. Zhu, Q.-f. Sun, T.-h. Lin, Andreev bound states and the $\pi$-junction transition in a superconductor/quantum-dot/superconductor system, Journal of Physics: Condensed Matter 13 (39) (2001) 8783.
  * [8] A. Kumar, M. Gaim, D. Steininger, A. L. Yeyati, A. Martín-Rodero, A. Hüttel, C. Strunk, Temperature dependence of andreev spectra in a superconducting carbon nanotube quantum dot, Physical Review B 89 (7) (2014) 075428\.
  * [9] J. Pillet, C. Quay, P. Morfin, C. Bena, A. L. Yeyati, P. Joyez, Andreev bound states in supercurrent-carrying carbon nanotubes revealed, Nature Physics 6 (12) (2010) 965–969.
  * [10] J. Kondo, Resistance minimum in dilute magnetic alloys, Progress of theoretical physics 32 (1) (1964) 37–49.
  * [11] B.-K. Kim, Y.-H. Ahn, J.-J. Kim, M.-S. Choi, M.-H. Bae, K. Kang, J. S. Lim, R. López, N. Kim, Transport measurement of andreev bound states in a kondo-correlated quantum dot, Physical review letters 110 (7) (2013) 076803.
  * [12] J. O. Island, R. Gaudenzi, J. De Bruijckere, E. Burzurí, C. Franco, M. Mas-Torrent, C. Rovira, J. Veciana, T. M. Klapwijk, R. Aguado, et al., Proximity-induced shiba states in a molecular junction, Physical Review Letters 118 (11) (2017) 117001.
  * [13] K. Franke, G. Schulze, J. Pascual, Competition of superconducting phenomena and kondo screening at the nanoscale, Science 332 (6032) (2011) 940–944.
  * [14] R. Žitko, M. Lee, R. López, R. Aguado, M.-S. Choi, Josephson current in strongly correlated double quantum dots, Phys. Rev. Lett. 105 (2010) 116803.  doi:10.1103/PhysRevLett.105.116803.
URL https://link.aps.org/doi/10.1103/PhysRevLett.105.116803

  * [15] J. A. Van Dam, Y. V. Nazarov, E. P. Bakkers, S. De Franceschi, L. P. Kouwenhoven, Supercurrent reversal in quantum dots, Nature 442 (7103) (2006) 667–670.
  * [16] S. De Franceschi, L. Kouwenhoven, C. Schönenberger, W. Wernsdorfer, Hybrid superconductor–quantum dot devices, Nature Nanotechnology 5 (10) (2010) 703.
  * [17] H. I. Jørgensen, T. Novotnỳ, K. Grove-Rasmussen, K. Flensberg, P. E. Lindelof, Critical current 0- $\pi$ transition in designed josephson quantum dot junctions, Nano letters 7 (8) (2007) 2441–2445.
  * [18] A. Y. Kasumov, K. Tsukagoshi, M. Kawamura, T. Kobayashi, Y. Aoyagi, K. Senba, T. Kodama, H. Nishikawa, I. Ikemoto, K. Kikuchi, et al., Proximity effect in a superconductor-metallofullerene-superconductor molecular junction, Physical Review B 72 (3) (2005) 033414.
  * [19] H. I. Jørgensen, T. Novotnỳ, K. Grove-Rasmussen, K. Flensberg, P. E. Lindelof, Critical current 0- $\pi$ transition in designed josephson quantum dot junctions, Nano letters 7 (8) (2007) 2441–2445.
  * [20] J.-D. Pillet, P. Joyez, M. Goffman, et al., Tunneling spectroscopy of a single quantum dot coupled to a superconductor: From kondo ridge to andreev bound states, Physical Review B 88 (4) (2013) 045101.
  * [21] M.-S. Choi, M. Lee, K. Kang, W. Belzig, Kondo effect and josephson current through a quantum dot between two superconductors, Physical Review B 70 (2) (2004) 020502.
  * [22] Y. Tanaka, A. Oguri, A. Hewson, Kondo effect in asymmetric josephson couplings through a quantum dot, New Journal of Physics 9 (5) (2007) 115.
  * [23] D. J. Luitz, F. F. Assaad, T. Novotnỳ, C. Karrasch, V. Meden, Understanding the josephson current through a kondo-correlated quantum dot, Physical review letters 108 (22) (2012) 227001.
  * [24] R. Delagrange, D. Luitz, R. Weil, A. Kasumov, V. Meden, H. Bouchiat, R. Deblock, Manipulating the magnetic state of a carbon nanotube josephson junction using the superconducting phase, Physical Review B 91 (24) (2015) 241401\.
  * [25] A. A. Clerk, V. Ambegaokar, Loss of $\pi$-junction behavior in an interacting impurity josephson junction, Physical Review B 61 (13) (2000) 9109.
  * [26] A. Martín-Rodero, A. Levy Yeyati, Josephson and andreev transport through quantum dots, Advances in Physics 60 (6) (2011) 899–958.
  * [27] S. Ishizaka, J. Sone, T. Ando, dc josephson current through a quantum dot coupled with superconducting leads, Physical Review B 52 (11) (1995) 8358.
  * [28] E. Vecino, A. Martín-Rodero, A. L. Yeyati, Josephson current through a correlated quantum level: Andreev states and $\pi$ junction behavior, Physical Review B 68 (3) (2003) 035105.
  * [29] J. F. Rentrop, S. G. Jakobs, V. Meden, Nonequilibrium transport through a josephson quantum dot, Physical Review B 89 (23) (2014) 235110.
  * [30] G. Rajput, R. Kumar, et al., Tunable josephson effect in hybrid parallel coupled double quantum dot-superconductor tunnel junction, Superlattices and Microstructures 73 (2014) 193–202.
  * [31] Y. Avishai, A. Golub, A. D. Zaikin, Superconductor-quantum dot-superconductor junction in the kondo regime, Physical Review B 67 (4) (2003) 041301.
  * [32] A. L. Yeyati, A. Martín-Rodero, E. Vecino, Nonequilibrium dynamics of andreev states in the kondo regime, Physical review letters 91 (26) (2003) 266802\.
  * [33] A. C. Hewson, The Kondo problem to heavy fermions, Vol. 2, Cambridge university press, 1997.


# Right-Handed Neutrino Dark Matter with Forbidden Annihilation

Yu Cheng<EMAIL_ADDRESS>Tsung-Dao Lee Institute & School of Physics and
Astronomy, Shanghai Jiao Tong University, China  Key Laboratory for Particle
Astrophysics and Cosmology (MOE) & Shanghai Key Laboratory for Particle
Physics and Cosmology, Shanghai Jiao Tong University, Shanghai 200240, China
Shao-Feng Ge<EMAIL_ADDRESS>Tsung-Dao Lee Institute & School of Physics and
Astronomy, Shanghai Jiao Tong University, China  Key Laboratory for Particle
Astrophysics and Cosmology (MOE) & Shanghai Key Laboratory for Particle
Physics and Cosmology, Shanghai Jiao Tong University, Shanghai 200240, China
Jie Sheng<EMAIL_ADDRESS>Tsung-Dao Lee Institute & School of Physics
and Astronomy, Shanghai Jiao Tong University, China  Key Laboratory for
Particle Astrophysics and Cosmology (MOE) & Shanghai Key Laboratory for
Particle Physics and Cosmology, Shanghai Jiao Tong University, Shanghai
200240, China  Tsutomu T. Yanagida<EMAIL_ADDRESS>Tsung-Dao Lee
Institute & School of Physics and Astronomy, Shanghai Jiao Tong University,
China  Key Laboratory for Particle Astrophysics and Cosmology (MOE) & Shanghai
Key Laboratory for Particle Physics and Cosmology, Shanghai Jiao Tong
University, Shanghai 200240, China

###### Abstract

The seesaw mechanism with three right-handed neutrinos has one as a well-
motivated dark matter candidate if stable and the other two can explain baryon
asymmetry via the thermal leptogenesis scenario. We explore the possibility of
introducing additional particles to make the right-handed neutrino dark matter
in thermal equilibrium and freeze out through a forbidden annihilation
channel. Nowadays in the Universe, this forbidden channel can be reactivated
by a strong gravitational potential such as the supermassive black hole in our
galaxy center. The Fermi-LAT gamma ray data and dark matter relic density
require this right-handed neutrino dark matter to have mass below $100\,$GeV
and the existence of an additional boson $\phi$ that can be tested at future
lepton colliders.

Introduction – The heavy right-handed Majorana neutrinos are widely considered
as a key ingredient beyond the Standard Model of particle physics. They
explain not only the observed tiny neutrino masses via the seesaw mechanism
Minkowski:1977sc ; Yanagida:1979as ; Yanagida:1979gs ; Gell-Mann:1979vob , but
also the Universe’s baryon asymmetry through leptogenesis Fukugita:1986hr .
Normally, we assume three heavy right-handed neutrinos (RHN). However, two
heavy Majorana neutrinos are sufficient to explain the observed neutrino
oscillation data and the baryon asymmetry in our Universe Frampton:2002qc ;
Raidal:2002xf ; Glashow:2003nk ; Barger:2003gt ; Ge:2010js . One remaining
right-handed neutrino can be either very heavy or very light Kusenko:2010ik ;
Drewes:2013gca ; Merle:2013gea ; Drewes:2016upu ; Dasgupta:2021ies . In this
letter, we explore the latter option that provides us with interesting
consequences at low energies.

If the remaining right-handed neutrino is as light as $\mathcal{O}$(100) GeV,
it can be a good candidate of dark matter (DM) Kusenko:2009up ; Merle:2013gea
; Drewes:2016upu ; Abazajian:2017tcc ; Boyarsky:2018tvu ; Dasgupta:2021ies .
However, we have strong constraints from cosmic microwave background (CMB) if
they annihilate into the SM particles and the subsequent electromagnetic
energy injection affects the ionization history of the universe Bean:2003kd ;
Chen:2003gz ; Hansen:2003yj ; Pierpaoli:2003rz ; Padmanabhan:2005es ;
Slatyer:2009yq ; Roszkowski:2017nbc ; Planck:2018nkj ; Cang:2020exa ;
Kawasaki:2021etm ; Liu:2023nct . Therefore, we consider a forbidden-type DM
DAgnolo:2015ujb ; Delgado:2016umt ; DAgnolo:2020mpt ; Wojcik:2021xki  whose
freeze-out production is through the annihilation channel to a heavier dark
sector partner or SM particles. In other words, DM annihilation is
kinematically forbidden. The only way out is the thermal energy that can
overcome the mass difference between the initial and final states. This can
happen at the early Universe but ceases when the temperature cools down with
the Universe expansion. Around the last scattering of CMB photons, the
temperature has dropped to eV scale which is negligibly small for any
reasonable forbidden DM model. So forbidden DM can avoid those CMB
constraints.

Without annihilation at present time, it is very difficult for observation to
verify the existence of this forbidden scenario. The indirect detection is
intrinsically forbidden. Interestingly, a strong gravitational source, such as
a supermassive Black Hole (SMBH), can reactivate the forbidden DM and make the
DM annihilate around it Cheng:2022esn . The subsequent decay of forbidden
channel final states into SM particles, especially visible photons, can give a
unique indirect detection signature today. The gamma ray from the DM forbidden
annihilation only appears in the vicinity around the SMBH, but not anywhere
else in the sky. A naive fit with the Fermi-LAT data for Sgr A* shows quite
good sensitivity on the thermally averaged annihilation cross section there.

Although the Fermi-LAT constraint is consistent with the DM relic density
requirement with naive assignments, a concrete model is necessary to
demonstrate the realistic possibility of constructing such a forbidden DM
theory. In this paper, we show our right-handed neutrino DM model and its
parameter space to generate the correct DM relic density and at the same time
provide reactivated annihilation signal. Comparing with the Fermi-LAT data, we
find the right-handed neutrino DM mass is below $100$ GeV and the presence of
an additional boson $\phi$ is required in the almost same mass region.
Finally, we point out that the possible test at future lepton collider.

The Right-Handed Neutrino Forbidden DM – We assume a $Z_{2}$ parity acting
only on the right-handed neutrino $N$ to decouple it from all the SM
particles. Then, $N$ becomes stable and a good DM candidate Cox:2017rgn  (If
it is the case, we have an anthropic reason for the presence of three families
of quarks and leptons Ibe:2016yfo .) However, its production is uncertain
since it decouples completely from the SM particles. It may be produced in the
early Universe through interactions with unknown heavy particles, but its
abundance is undetermined by low energy physics.

A possible solution is to introduce a fermion $\chi$ and a boson $\phi$ to
couple with the DM $N$, $\phi N\chi$. Notice here that all fermions, $N$ and
$\chi$, are left-handed two component Weyl fermions. The DM neutrino $N$ is
assumed to carry odd parity under the $Z_{2}$. If $\phi$ is even (and $\chi$
is odd), it can couple to a pair of the SM Higgs bosons, $H$ and
$H^{\dagger}$, via $\phi H^{\dagger}H$. Thus, the interaction Lagrangian is,

$\mathcal{L}_{\rm int}=(y\phi N\chi+h.c.)+\lambda m_{\phi}\phi H^{\dagger}H,$
(1)

where $m_{\phi}$ is the mass of $\phi$. The coupling $\lambda$ with the SM
Higgs boson renders $\phi$ in thermal equilibrium with the SM particles in the
early Universe.

To maintain equilibrium, the decay rate of $\phi$ to SM fermions,
$\Gamma\sim\frac{1}{32\pi}\lambda^{2}m^{2}_{f}m^{3}_{\phi}/(m_{h}^{2}-m_{\phi}^{2})^{2}$,
should be larger than the Hubble constant $H\propto T_{f}^{2}/M_{\rm pl}$
where $M_{\rm Pl}$ is the Planck mass. Suppose the $\phi$ decoupling happens
around $T_{f}\sim m_{\phi}/25$, thermal equilibrium requires
$\lambda\gtrsim(10^{-8}\sim 10^{-7})$ depending on the scalar mass $m_{\phi}$.
Even lower limit is possible if $m_{\phi}$ approaches the Higgs mass to reduce
the denominator $m^{2}_{\phi}-m^{2}_{h}$ of the decay width formula. Further
through the Yukawa term and the resulting $NN\leftrightarrow\phi\phi$
scatterings, $N$ can also get in the thermal equilibrium.

All these new particles are singlets of the SM gauge group and have mass
terms,

$\mathcal{L}=\frac{m_{N}}{2}NN+\frac{m_{\chi}}{2}\chi\chi+\frac{m_{\phi}^{2}}{2}\phi^{2}.$
(2)

For a forbidden scenario, the DM mass should be lighter than its annihilation
final-state particles, $m_{N}<m_{\phi}$. Otherwise, the RHN $N$ can pairly
annihilate into a pair of $\phi$ at the late time and then $\phi$ decays to SM
particles through the mixing with the SM Higgs boson. For $\phi$ mass larger
than $\mathcal{O}(1)$ GeV, the dominated decay channel is $\phi\rightarrow
f\bar{f}$, where $f=b,\tau$. If $\phi$ mass is $\lesssim\mathcal{O}(1)$ GeV,
$\phi$ can still decay through the di-photon channel
$\phi\rightarrow\gamma\gamma$ that is induced by triangle diagrams in the
similar way as the SM Higgs boson Workman:2022ynf . These processes are
strongly constrained by the CMB Elor:2015bho  and high energy gamma ray data
Fortes:2022cxp . Therefore, we concentrate our discussion on the case that the
boson $\phi$ is heavier than the DM $N$, $m_{\phi}>m_{N}$. Although the DM $N$
may annihilate into SM fermions by exchanging a virtual Higgs,
$NN\rightarrow\phi+h^{*}$ and $h^{*}\rightarrow f\bar{f}$, this three-body
final state is naturally suppressed by a phase factor $1/16\pi^{2}\sim
10^{-2}$ and the Yukawa coupling between Higgs and fermions $y^{2}\lesssim
10^{-3}$. In total, it has a $\mathcal{O}(10^{-5})$ suppression relative to
the two body annihilation so that the CMB constraint can hardly apply. We also
consider that the mass of $\chi$ is sufficiently heavy,
$m_{\chi}>m_{\phi}+m_{N}$, such that $\chi$ can decay and does not contribute
to the DM relic density. Below we focus on the density evolution of the right-
handed neutrino forbidden DM $N$.

In the beginning, the DM $N$ is also in thermal equilibrium. It gradually
freezes out through the forbidden annihilation channel
$N+N\leftrightarrow\phi+\phi$ by exchanging a $t$-channel $\chi$. The number
density evolution is governed by the Boltzmann equation,

$\dot{n}_{N}+3Hn_{N}=-\left\langle\sigma_{NN}v\right\rangle\left[n_{N}^{2}-\left(n^{\rm
eq}_{N}\right)^{2}\right],$ (3)

where $n_{N}$ ($n_{N}^{\rm eq})$ is the number density of right-handed
neutrino $N$ (in equilibrium). Since the kinematic threshold of the forbidden
channel cuts the allowed range of the center-of-mass energy $\sqrt{s}$, the
total cross section cannot easily find an analytical formula. It is more
convenient to evaluate the cross section $\langle\sigma_{\phi\phi}v\rangle$ of
the reversed process $\phi\phi\rightarrow NN$ first and express
$\left\langle\sigma_{NN}v\right\rangle$ according to detailed balance
DAgnolo:2015ujb ,

$\left\langle\sigma_{NN}v\right\rangle\equiv\left(\frac{n_{\phi}^{\rm
eq}}{n_{N}^{\rm
eq}}\right)^{2}\left\langle\sigma_{\phi\phi}v\right\rangle=\frac{(1+\delta)^{3}}{4}e^{-2x\delta}\left\langle\sigma_{\phi\phi}v\right\rangle.$
(4)

Being a function of the mass difference $\delta\equiv(m_{\phi}-m_{N})/m_{N}$
and temperature parameter $x\equiv m_{N}/T$, the exponential suppression
$e^{-2x\delta}$ comes from the ratio of number densities in equilibrium,
$n_{\rm eq}=g(mT/2\pi)^{3/2}\exp(-m/T)$, where $g=2$ and 1 for $N$ and $\phi$,
respectively. The Boltzmann equation can be solved semi-analytically to obtain
the rescaled number density $Y_{N}\equiv n_{N}/T^{3}$,

$Y_{N}=\frac{\left[(a-b\delta)\frac{1}{x_{f}}e^{-2\delta
x_{f}}g_{\delta}\left(x_{f}\right)+\frac{b}{2x_{f}^{2}}e^{-2\delta
x_{f}}\right]^{-1}}{\sqrt{\frac{\pi}{45}}m_{N}M_{\mathrm{Pl}}(1+\delta)^{3}g^{s}_{*}/\sqrt{g_{*}}}$
(5)

where $x_{f}$ indicates the freeze-out temperature and $g^{s}_{*}$ ($g_{*}$)
the entropy (total) degrees of freedom then. For convenience, we have defined
$g_{\delta}(x_{f})$ $\equiv$ $1-2\delta x_{f}$ $e^{2\delta x_{f}}$
$\int_{2\delta x_{f}}^{\infty}$ $t^{-1}$ $e^{-t}$ $dt$. The parameters $a$ and
$b$ are the expansion coefficients of the thermally averaged cross section
$\left\langle\sigma_{\phi\phi}v\right\rangle$ for the $s$ and $p$ waves,
$\left\langle\sigma_{\phi\phi}v\right\rangle\equiv a+bx^{-1}$ Gondolo:1990dk .
The DM relic density is related to $Y_{N}$ as,

$\rho_{\mathrm{DM}}=m_{N}s_{0}Y_{N},\quad\Omega_{\mathrm{DM}}h^{2}=\frac{\rho_{\mathrm{DM}}}{\rho_{c}/h^{2}},$
(6)

where $s_{0}=2891.2$ cm-3 is today’s entropy density and $\rho_{c}=1.05\times
10^{-5}h^{2}$ GeV/cm3 is the critical density of the Universe.

We also solve the Boltzmann equation numerically with MicrOMEGAs
Belanger:2020gnr . The allowed parameter space to produce the correct relic
density $\Omega_{\rm DM}h^{2}=0.12$ is shown in Fig. 1. The two key mass
parameters of the forbidden channel are shown in terms of the right-handed
neutrino mass $m_{N}$ and the relative mass difference $\delta$ as horizontal
and vertical axes, respectively. The whole $m_{N}$–$\delta$ parameter space is
all allowed by tuning the Yukawa coupling $y$ to match the DM relic density.
In other words, the value of $y$ is given as a function of $m_{N}$, $\delta$,
and $\Omega_{\rm DM}$. Although the partner fermion $\chi$ also has an
adjustable mass $m_{\chi}$ that can enter the annihilation cross section,
tuning its value affects the DM relic density in the early Universe and the
reactivated forbidden annihilation at current days simultaneously. This
resembles the effect of the Yukawa coupling $y$. Tuning the value of either
Yukawa coupling $y$ or $m_{\chi}$ is sufficient for our purpose. For
simplicity, we take fixed $m_{\chi}=2(m_{\phi}+m_{N})$ for illustration. Fig.
1 indicates that the DM relic density naturally limits its mass within the
range of $m_{N}\lesssim 100$ GeV. A smaller coupling constant $y$ prefers
smaller DM mass $m_{N}$ and mass difference $\delta$. This is because when the
annihilation cross section becomes smaller, DM freezes out earlier and its
number density becomes larger. Since the DM relic density is inversely
proportional to the annihilation cross section $\braket{\sigma v}$, it scales
with $m_{N}^{2}/y^{4}$. A smaller coupling requires a lighter DM. On the other
hand, a smaller mass difference can increase the annihilation cross section to
compensate the effect of a small coupling.

The Forbidden DM Signature – The forbidden DM scenario has an intrinsic
difficulty of being probed via indirect detection. Nowadays, the DM Boltzmann
distribution with typical $\mathcal{O}(100)$ km/s velocity in our galaxy can
no longer support the forbidden channel. Such a non-relativistic velocity
$v\sim 10^{-3}$ can only overcome the mass threshold $\delta\sim v^{2}\sim
10^{-6}$. Although direct detection and collider experiments may provide
complementary searches and even measure the DM mass, the forbidden nature of
DM is still missing. In other words, the theoretical models of forbidden DM
can in principle exist but is not testable. This is a very serious issue.

To reopen the forbidden channel, it is necessary to accelerate the DM
particles. Around a black hole, especially the SMBH in the galaxy center, a
particle can be accelerated to become relativistic because of the strong
gravitational force Cheng:2022esn . This can help the DM $N$ to annihilate
into heavier dark sector particle $\phi$ near the SMBH. The unstable $\phi$
then decays to SM particles through its mixing with the Higgs particle. The
decay products are either gamma photons or charged particles. While the di-
photon decay channel gives mono-energetic gamma rays in the $\phi$ rest frame,
the charged particles produce photons through final-state radiation or meson
decays with continuous spectrum. These allow using gamma ray that can point
straightly back to the SMBH as signal for testing our model of right-handed
neutrino DM.

Figure 1:  The allowed parameter space contour to obtain the DM relic density
through freeze out. After fixing non-DM fermion mass
$m_{\chi}=2(m_{N}+m_{\phi})$, there are still three model parameters, the
right-handed neutrino DM mass $m_{N}$ (horizontal axis), the mass difference
$\delta$ (vertical axis), and the Yukawa coupling constant $y$ (contours
according to the color bar). Red lines show the 95% limit provided by the
FermiLAT gamma ray data. Both cases with fixed (solid red) and unfixed (dashed
red) background are shown for comparison. The almost vertical black dashed
line indicates the maximial $m_{N}$ that allows Higgs invisible decay
$h\rightarrow N+\chi$.

The gamma ray intensity highly depends on the DM density and velocity
dispersion profiles around the SMBH. In our scenario, DM $N$ has no self-
interaction. Thus, its density follows the CDM density profile Fields:2014pia
. For the innermost region, $r<4GM\equiv r_{0}$ where $M$ is the BH mass, all
particles are attracted to fall into the SMBH so that the DM density is $0$
therein. On the other hand, the DM halo simply follows the NFW profile when
the gravitational influence of the SMBH no longer dominates which roughy
happens at $r_{b}\equiv 0.2GM/v_{0}^{2}$ with $v_{0}$ being the DM velocity
dispersion there. In between, a DM density spike forms,

$\rho(r)=\begin{cases}0,&r<r_{0},\text{ (Capture Region),}\\\
\frac{\rho_{\mathrm{sp}}(r)\rho_{\mathrm{in}}(t,r)}{\rho_{\mathrm{sp}}(r)+\rho_{\mathrm{in}}(t,r)},&r_{0}\leq
r<r_{b},\text{ (Spike), }\\\ \text{ NFW Profile},&r>r_{b},\text{ (Halo).
}\end{cases}$ (7)

Naively thinking, the DM density profiles keeps increasing when going towards
the SMBH. However, the DM density cannot increase forever. In addition to
self-interaction, annihilation can also suppresses the growth of DM density.
This naturally puts an upper bound on the DM density profile and forms a
density plateau around the SMBH. Going inward from the spike region boundary
$r_{b}$, the DM density profile increases very fast initially as $\rho_{\rm
sp}(r)\equiv\rho_{b}(r_{b}/r)^{\gamma_{\rm sp}}$ with $\gamma_{\rm
sp}\equiv(9-2\gamma_{c})/(4-\gamma_{c})$. The density coefficient
$\rho_{b}\equiv\rho_{D}(D/r_{b})^{\gamma_{c}}$ with $\gamma_{c}=1$ is scaled
according to the NFW profile Navarro:1995iw , $\rho(r)\propto 1/r$ for $r\ll
26$ kpc Alvarez:2020fyo ; Abazajian:2020tww , from the density around the
Solar system with $\rho_{D}=0.3\,$ GeV cm-3 and $D=8.5$ kpc. When reaching the
SMBH vicinity, $r\gtrsim r_{0}$, the DM density increases much slower as
$\rho_{\rm in}(r,t)\equiv\rho_{\rm ann}(t)(r/r_{\rm in})^{-\gamma_{\rm in}}$
with $\gamma_{\rm in}=1/2$. The density $\rho_{\rm ann}$ is the so-called
annihilation plateau density $\rho_{\rm ann}\equiv
m_{N}/\braket{\sigma_{NN}v}t$. With $\rho=\rho_{\rm sp}\rho_{\rm
in}/(\rho_{\rm sp}+\rho_{\rm in})$, its value approaches $\rho_{\rm sp}$ for
$\rho_{\rm sp}\ll\rho_{\rm in}$ and $\rho_{\rm in}$ in the opposite. Since
$\rho_{\rm sp}(r)\propto r^{-7/3}$ and $\rho_{\rm in}(r)\propto r^{-1/2}$, the
DM density profile almost becomes a plateau in the inner spike region. The
division $r_{\rm in}$ at $\rho_{\rm sp}=\rho_{\rm ann}$ is determined by the
DM annihilation cross section $\langle\sigma_{NN}v\rangle$.

In its vicinity, the SMBH dominates the gravitational potential and
consequently the DM velocity dispersion follows a simple scaling
$v_{d}^{2}(r)\sim GM/r$. In reality, the DM particles inside the spike region
are thermalized to a Juttner distribution Juttner:1911 ; DeGroot:1980dk
$P_{r}(V_{r},V_{c},x(r))$ where $V_{r}$ is the relative velocity and $V_{c}$
is the centre-of-mass velocity of a two-particle collision system
Cheng:2022esn . To get the energy spectrum of the gamma production rate per
volume $d\Phi_{\gamma}(r)/dE_{\gamma}$, one needs to integrate over the DM
distribution from the threshold relative velocity $V^{\rm
th}_{r}\equiv[1-1/(1-2m_{F}^{2}/m_{\chi}^{2})^{2}]^{1/2}$,

$\frac{dF_{\gamma}}{dE_{\gamma}}(r)=\int_{V^{\rm
th}}^{1}dV_{r}P_{r}\left(V_{r},V_{c},x(r)\right)\sigma
V_{r}\frac{dN_{\gamma}}{dE_{\gamma}}(V_{r},V_{c}),$ (8)

where $x(r)$ has radius dependence through the DM temperature
$T=\frac{1}{2}m_{N}v_{d}^{2}(r)$ and $dN_{\gamma}/dE_{\gamma}$ is the boosted
photon spectrum from $\phi$ decay. For the mass range $m_{\phi}>10$ GeV
considered in this paper, the dominate channels are $\phi\rightarrow
b\bar{b}$, $\tau\bar{\tau}$, and $gg$. We use the PPPC4DMID package
Cirelli:2010xx  to generate the photon spectra of $\phi$ decay at rest and
HDECAY Djouadi:1997yw  to obtain the branching ratios. After that, we boost
the spectrum to the galaxy frame Elor:2015bho ; Elor:2015tva .

The observed gamma ray flux is an integrated total result over the radius from
$r_{0}=4GM$ to some upper limit $r_{B}$,

$\frac{d\Phi_{\gamma}}{dE_{\gamma}}=\frac{1}{4\pi
D^{2}}\frac{1}{2m_{\chi}^{2}}\int_{4GM}^{r_{B}}4\pi
r^{2}dr\rho^{2}(r)\frac{dF_{\gamma}}{dE_{\gamma}}(r).$ (9)

In principle, the integration upper limit $r_{B}$ should be as large as
possible to include all contributions. However, those DM particles at larger
radius have a smaller velocity dispersion and consequently only the tail of
the Juttner distribution contributes. To make efficient numerical evaluation,
we only include those regions that,

$\int^{1}_{V^{\rm
th}_{r}}dV_{r}\int_{0}^{1}dV_{c}P_{r}(V_{r},V_{c},x(r_{B}))\equiv 1\%.$ (10)

This means that only less than 1% of the DM phase space can contribute at
radius larger than $r_{B}$.

Indirect Detection with Fermi-LAT – The photon flux from the forbidden
annihilation is localized around the SMBH. Uniquely identifying the right-
handed neutrino DM scenario requires not just observing the photon spectrum
but also good angular resolution. The Fermi-LAT satellite observatory is an
all-sky $\gamma$-ray telescope with excellent energy and angular resolutions
on the $\gamma$-ray point source from the GC Fermi-LAT:2009ihh ; Fermi-
LAT:2016uux . We analyze a square region of $10^{\circ}\times 10^{\circ}$
around Sgr $A^{*}$ and find several $\gamma$-ray point sources, of which 4FGL
J1745.6-2859 is the brightest and closest one to Sgr $A^{*}$ in the Fourth
catalog of Fermi-LAT sources (4FGL) Fermi-LAT:2019yla ; Ballet:2020hze . This
point source is usually considered as the manifestation of Sgr $A^{*}$ in the
MeV-to-GeV range Cafardo:2021pqs . We use 14 years of Fermi-LAT data from
August 4, 2008 to October 26, 2022 in this paper. To be specific, the Pass 8
SOURCE-class events from 100 MeV to 1000 GeV are binned to a pixel size of
$0.08^{\circ}$.

A universal model can describe the $\gamma$-ray spectrum from point sources.
To be specific, this universal 4FGL J1745.6-2859 spectral model is a log-
parabola in the 4FGL Catalog Fermi-LAT:2019yla ; Ballet:2020hze ,

$\frac{dN}{dE}=N_{0}\left(\frac{E}{E_{0}}\right)^{-\alpha-\beta\log(E/E_{0})},$
(11)

where $N_{0}$ is normalization, $E_{0}$ a scale parameter, $\alpha$ the
spectral slope at $E_{0}$, and $\beta$ the curvature of the spectrum. Since
$E_{0}$ does not vary much, its value is fixed to 4074 MeV in our fit while
the other three can freely adjust.

For a given DM mass $m_{N}$ and mass difference $\delta$, the Yukawa coupling
constant $y$ is fixed by the DM relic density as explained earlier. This means
that the predicted photon flux around the SMBH is also determined for every
parameter space point in Fig. 1. We fit the Fermi-LAT data with the following
$\chi^{2}$ function,

$\chi^{2}\equiv\sum_{i}^{\text{bins
}}\left(\frac{N_{i}^{\mathrm{data}}-N_{i}^{\mathrm{bkg}}-N_{i}^{\mathrm{NP}}}{\sigma_{i}}\right)^{2}.$
(12)

In each bin, $N_{i}^{\mathrm{data}}$ is the data point, $N_{i}^{\mathrm{bkg}}$
the background, $N_{i}^{\mathrm{NP}}$ the expected events due to new physics
(NP) contribution, and $\sigma_{i}$ the error bar for each data point.

Fig. 1 shows the 95% C.L. limit as red lines and the region above are
excluded. Since the realistic background is not determined and the parameters
in Eq. (11) are free, we do the fit in two possible ways. We first fix the
background to its best-fit values obtained from the background-only hypothesis
with $N_{0}=2.519\times 10^{-12}$, $\alpha=2.489$, and $\beta=0.1553$. In this
case, the background model can fit data very well and consequently the
resulting limit on the signal parameters (red solid line) are quite stringent.
The allowed DM mass is below $\lesssim 60$ GeV and the mass difference
$\delta$ cannot be larger than $\mathcal{O}(10^{-5})$. For comparison, we also
fit the data by changing both the prediction and background parameters. The
corresponding result (red dashed line) allows much larger parameter space with
a mass difference as large as $10^{-3}$. Both methods require the mass
difference within a quite narrow range at the leave of $\mathcal{O}(10^{-5})$
and $\mathcal{O}(10^{-3})$, respectively for DM mass $m_{N}\gtrsim 30$ GeV.

Conclusion and Discussion – In this paper, we propose a model with dark sector
including two Majorana fermions $N$ and $\chi$, as well as a Higgs-portal
scalar particle $\phi$. Among them, the right-handed neutrino $N$ is the
lightest one and hence can serve as DM. Its freeze-out production in the early
Universe is through the forbidden annihilation channel
$N+N\leftrightarrow\phi+\phi$. Although this annihilation channel is nowadays
forbidden, it can reopen around a strong gravitation source. Thus, the SMBH
located at the galaxy center can accelerate $N$ to make it annihilate into its
heavier partner $\phi$. The subsequent decay of $\phi$ to SM particles can
bring detectable photon flux. The current Fermi-LAT data requires the right-
handed neutrino mass $m_{N}$ smaller than $100$ GeV, and the mass difference
$\mathcal{O}(10^{-5})\lesssim\delta\lesssim\mathcal{O}(10^{-3})$.

Such a light right-handed neutrino can be tested at future lepton colliders
such as CEPC CEPCStudyGroup:2018ghi , ILC ILC:2013jhg , FCC-ee FCC:2018evy ,
and CLIC Linssen:2012hp . If the masses of $N$ and $\chi$ satisfy
$m_{\chi}+m_{N}<m_{h}$, the SM Higgs particle can decay invisibly,
$h\rightarrow\chi+N$. For illustration, the projected sensitivity
$Br_{\mathrm{inv}}^{\mathrm{BSM}}<0.3\%$ at CEPC translates to
$\lambda\lesssim 5\times 10^{-2}$. With our parameter assignments,
$m_{\phi}>m_{N}$ and $m_{\chi}>2(m_{N}+m_{\phi})$, this Higgs invisible decay
constraint applies for $m_{N}<25$ GeV as shown with verticle black dashed line
in Fig. 1. Besides, the dark sector particle $\phi$ can also be directly
produced at future Higgs factories through Higgsstrahlung-like process
$e^{+}+e^{-}\rightarrow Z+\phi$ Ko:2016xwd ; Heng:2017jgh ; Liu:2017lpo ;
Kamon:2017yfx ; Azevedo:2018oxv ; Grzadkowski:2020frj ; Haghighat:2022qyh .

In addition to the current model assignments discussed in this paper, there is
another possibility of removing the dark sector particle $\chi$ and assigning
the right-handed neutrino DM $N$ play its role. Then the forbidden
annihilation $NN\rightarrow\phi\phi$ arises from a much simpler flavor-
conserving Yukawa coupling $y\phi NN$. One side effect of this model is
introducing DM self-scattering $NN\leftrightarrow NN$ by exchanging a scalar
mediator $\phi$. Consequently, DM has much a smaller density profile in the
spike region. This leads to a weaker constraint on the annihilation cross
section FuturePaper .

## Acknowledgements

SFG is supported by the Double First Class start-up fund (WF220442604) and
National Natural Science Foundation of China (No. 12090064), and Chinese
Academy of Sciences Center for Excellence in Particle Physics (CCEPP). T. T.
Y. is supported by the China Grant for Talent Scientific Start-Up Project and
by Natural Science Foundation of China (NSFC) under grant No. 12175134, JSPS
Grant-in-Aid for Scientific Research Grants No. 19H05810, and World Premier
International Research Center Initiative (WPI Initiative), MEXT, Japan. Both
SFG and T. T. Y. are affiliated members of Kavli IPMU, University of Tokyo.

## References

  * (1) P. Minkowski, “$\mu\to e\gamma$ at a Rate of One Out of $10^{9}$ Muon Decays?,” Phys. Lett. B 67 (1977) 421–428.
  * (2) T. Yanagida, “Horizontal gauge symmetry and masses of neutrinos,” Conf. Proc. C 7902131 (1979) 95–99.
  * (3) T. Yanagida, “Horizontal Symmetry and Mass of the Top Quark,” Phys. Rev. D 20 (1979) 2986.
  * (4) M. Gell-Mann, P. Ramond, and R. Slansky, “Complex Spinors and Unified Theories,” Conf. Proc. C 790927 (1979) 315–321, [arXiv:1306.4669 [hep-th]].
  * (5) M. Fukugita and T. Yanagida, “Baryogenesis Without Grand Unification,” Phys. Lett. B 174 (1986) 45–47.
  * (6) P. H. Frampton, S. L. Glashow, and T. Yanagida, “Cosmological sign of neutrino CP violation,” Phys. Lett. B 548 (2002) 119–121, [arXiv:hep-ph/0208157].
  * (7) M. Raidal and A. Strumia, “Predictions of the most minimal seesaw model,” Phys. Lett. B 553 (2003) 72–78, [arXiv:hep-ph/0210021].
  * (8) S. L. Glashow, “Fact and fancy in neutrino physics. 2.,” in 10th International Workshop on Neutrino Telescopes, pp. 611–617.  6, 2003.  [arXiv:hep-ph/0306100].
  * (9) V. Barger, D. A. Dicus, H.-J. He, and T.-j. Li, “Structure of cosmological CP violation via neutrino seesaw,” Phys. Lett. B 583 (2004) 173–185, [arXiv:hep-ph/0310278].
  * (10) S.-F. Ge, H.-J. He, and F.-R. Yin, “Common Origin of Soft mu-tau and CP Breaking in Neutrino Seesaw and the Origin of Matter,” JCAP 05 (2010) 017, [arXiv:1001.0940 [hep-ph]].
  * (11) A. Kusenko, F. Takahashi, and T. T. Yanagida, “Dark Matter from Split Seesaw,” Phys. Lett. B 693 (2010) 144–148, [arXiv:1006.1731 [hep-ph]].
  * (12) M. Drewes, “The Phenomenology of Right Handed Neutrinos,” Int. J. Mod. Phys. E 22 (2013) 1330019, [arXiv:1303.6912 [hep-ph]].
  * (13) A. Merle, “keV Neutrino Model Building,” Int. J. Mod. Phys. D 22 (2013) 1330020, [arXiv:1302.2625 [hep-ph]].
  * (14) M. Drewes et al., “A White Paper on keV Sterile Neutrino Dark Matter,” JCAP 01 (2017) 025, [arXiv:1602.04816 [hep-ph]].
  * (15) B. Dasgupta and J. Kopp, “Sterile Neutrinos,” Phys. Rept. 928 (2021) 1–63, [arXiv:2106.05913 [hep-ph]].
  * (16) A. Kusenko, “Sterile neutrinos: The Dark side of the light fermions,” Phys. Rept. 481 (2009) 1–28, [arXiv:0906.2968 [hep-ph]].
  * (17) K. N. Abazajian, “Sterile neutrinos in cosmology,” Phys. Rept. 711-712 (2017) 1–28, [arXiv:1705.01837 [hep-ph]].
  * (18) A. Boyarsky, M. Drewes, T. Lasserre, S. Mertens, and O. Ruchayskiy, “Sterile neutrino Dark Matter,” Prog. Part. Nucl. Phys. 104 (2019) 1–45, [arXiv:1807.07938 [hep-ph]].
  * (19) R. Bean, A. Melchiorri, and J. Silk, “Recombining WMAP: Constraints on ionizing and resonance radiation at recombination,” Phys. Rev. D 68 (2003) 083501, [arXiv:astro-ph/0306357].
  * (20) X.-L. Chen and M. Kamionkowski, “Particle decays during the cosmic dark ages,” Phys. Rev. D 70 (2004) 043502, [arXiv:astro-ph/0310473].
  * (21) S. H. Hansen and Z. Haiman, “Do we need stars to reionize the universe at high redshifts? Early reionization by decaying heavy sterile neutrinos,” Astrophys. J. 600 (2004) 26–31, [arXiv:astro-ph/0305126].
  * (22) E. Pierpaoli, “Decaying particles and the reionization history of the universe,” Phys. Rev. Lett. 92 (2004) 031301, [arXiv:astro-ph/0310375].
  * (23) N. Padmanabhan and D. P. Finkbeiner, “Detecting dark matter annihilation with CMB polarization: Signatures and experimental prospects,” Phys. Rev. D 72 (2005) 023508, [arXiv:astro-ph/0503486].
  * (24) T. R. Slatyer, N. Padmanabhan, and D. P. Finkbeiner, “CMB Constraints on WIMP Annihilation: Energy Absorption During the Recombination Epoch,” Phys. Rev. D 80 (2009) 043526, [arXiv:0906.1197 [astro-ph.CO]].
  * (25) L. Roszkowski, E. M. Sessolo, and S. Trojanowski, “WIMP dark matter candidates and searches—current status and future prospects,” Rept. Prog. Phys. 81 no. 6, (2018) 066201, [arXiv:1707.06277 [hep-ph]].
  * (26) Planck Collaboration, N. Aghanim et al., “Planck 2018 results. I. Overview and the cosmological legacy of Planck,” Astron. Astrophys. 641 (2020) A1, [arXiv:1807.06205 [astro-ph.CO]].
  * (27) J. Cang, Y. Gao, and Y.-Z. Ma, “Probing dark matter with future CMB measurements,” Phys. Rev. D 102 no. 10, (2020) 103005, [arXiv:2002.03380 [astro-ph.CO]].
  * (28) M. Kawasaki, H. Nakatsuka, K. Nakayama, and T. Sekiguchi, “Revisiting CMB constraints on dark matter annihilation,” JCAP 12 no. 12, (2021) 015, [arXiv:2105.08334 [astro-ph.CO]].
  * (29) H. Liu, W. Qin, G. W. Ridgway, and T. R. Slatyer, “Exotic energy injection in the early universe II: CMB spectral distortions and constraints on light dark matter,” [arXiv:2303.07370 [astro-ph.CO]].
  * (30) R. T. D’Agnolo and J. T. Ruderman, “Light Dark Matter from Forbidden Channels,” Phys. Rev. Lett. 115 no. 6, (2015) 061301, [arXiv:1505.07107 [hep-ph]].
  * (31) A. Delgado, A. Martin, and N. Raj, “Forbidden Dark Matter at the Weak Scale via the Top Portal,” Phys. Rev. D 95 no. 3, (2017) 035002, [arXiv:1608.05345 [hep-ph]].
  * (32) R. T. D’Agnolo, D. Liu, J. T. Ruderman, and P.-J. Wang, “Forbidden dark matter annihilations into Standard Model particles,” JHEP 06 (2021) 103, [arXiv:2012.11766 [hep-ph]].
  * (33) G. N. Wojcik and T. G. Rizzo, “Forbidden scalar dark matter and dark Higgses,” JHEP 04 (2022) 033, [arXiv:2109.07369 [hep-ph]].
  * (34) Y. Cheng, S.-F. Ge, X.-G. He, and J. Sheng, “Forbidden Dark Matter Combusted Around Supermassive Black Hole,” [arXiv:2211.05643 [hep-ph]].
  * (35) P. Cox, C. Han, and T. T. Yanagida, “Right-handed Neutrino Dark Matter in a U(1) Extension of the Standard Model,” JCAP 01 (2018) 029, [arXiv:1710.01585 [hep-ph]].
  * (36) M. Ibe, A. Kusenko, and T. T. Yanagida, “Why three generations?,” Phys. Lett. B 758 (2016) 365–369, [arXiv:1602.03003 [hep-ph]].
  * (37) Particle Data Group Collaboration, R. L. Workman et al., “Review of Particle Physics,” PTEP 2022 (2022) 083C01.
  * (38) G. Elor, N. L. Rodd, T. R. Slatyer, and W. Xue, “Model-Independent Indirect Detection Constraints on Hidden Sector Dark Matter,” JCAP 06 (2016) 024, [arXiv:1511.08787 [hep-ph]].
  * (39) G. Fortes, F. Queiroz, C. Siqueira, and A. Viana, “Present and Future constraints on Secluded Dark Matter in the Galactic Halo with TeV Gamma-ray Observatories,” [arXiv:2212.05075 [hep-ph]].
  * (40) P. Gondolo and G. Gelmini, “Cosmic abundances of stable particles: Improved analysis,” Nucl. Phys. B 360 (1991) 145–179.
  * (41) G. Belanger, A. Mjallal, and A. Pukhov, “Recasting direct detection limits within micrOMEGAs and implication for non-standard Dark Matter scenarios,” Eur. Phys. J. C 81 no. 3, (2021) 239, [arXiv:2003.08621 [hep-ph]].
  * (42) B. D. Fields, S. L. Shapiro, and J. Shelton, “Galactic Center Gamma-Ray Excess from Dark Matter Annihilation: Is There A Black Hole Spike?,” Phys. Rev. Lett. 113 (2014) 151302, [arXiv:1406.4856 [astro-ph.HE]].
  * (43) J. F. Navarro, C. S. Frenk, and S. D. M. White, “The Structure of cold dark matter halos,” Astrophys. J. 462 (1996) 563–575, [arXiv:astro-ph/9508025].
  * (44) G. Alvarez and H.-B. Yu, “Density spikes near black holes in self-interacting dark matter halos and indirect detection constraints,” Phys. Rev. D 104 no. 4, (2021) 043013, [arXiv:2012.15050 [hep-ph]].
  * (45) K. N. Abazajian, S. Horiuchi, M. Kaplinghat, R. E. Keeley, and O. Macias, “Strong constraints on thermal relic dark matter from Fermi-LAT observations of the Galactic Center,” Phys. Rev. D 102 no. 4, (2020) 043012, [arXiv:2003.10416 [hep-ph]].
  * (46) F. Jüttner., “Das Maxwellsche Gesetz der Geschwindigkeitsverteilung in der Relativtheorie,” Ann. Phys. 339 (1911) 856.
  * (47) S. R. De Groot, Relativistic Kinetic Theory. Principles and Applications.  1980\.
  * (48) M. Cirelli, G. Corcella, A. Hektor, G. Hutsi, M. Kadastik, P. Panci, M. Raidal, F. Sala, and A. Strumia, “PPPC 4 DM ID: A Poor Particle Physicist Cookbook for Dark Matter Indirect Detection,” JCAP 03 (2011) 051, [arXiv:1012.4515 [hep-ph]]. [Erratum: JCAP 10, E01 (2012)].
  * (49) A. Djouadi, J. Kalinowski, and M. Spira, “HDECAY: A Program for Higgs boson decays in the standard model and its supersymmetric extension,” Comput. Phys. Commun. 108 (1998) 56–74, [arXiv:hep-ph/9704448].
  * (50) G. Elor, N. L. Rodd, and T. R. Slatyer, “Multistep cascade annihilations of dark matter and the Galactic Center excess,” Phys. Rev. D 91 (2015) 103531, [arXiv:1503.01773 [hep-ph]].
  * (51) Fermi-LAT Collaboration, W. B. Atwood et al., “The Large Area Telescope on the Fermi Gamma-ray Space Telescope Mission,” Astrophys. J. 697 (2009) 1071–1102, [arXiv:0902.1089 [astro-ph.IM]].
  * (52) Fermi-LAT, DES Collaboration, A. Albert et al., “Searching for Dark Matter Annihilation in Recently Discovered Milky Way Satellites with Fermi-LAT,” Astrophys. J. 834 no. 2, (2017) 110, [arXiv:1611.03184 [astro-ph.HE]].
  * (53) Fermi-LAT Collaboration, S. Abdollahi et al., “$Fermi$ Large Area Telescope Fourth Source Catalog,” Astrophys. J. Suppl. 247 no. 1, (2020) 33, [arXiv:1902.10045 [astro-ph.HE]].
  * (54) Fermi-LAT Collaboration, J. Ballet, T. H. Burnett, S. W. Digel, and B. Lott, “Fermi Large Area Telescope Fourth Source Catalog Data Release 2,” [arXiv:2005.11208 [astro-ph.HE]].
  * (55) Fermi-LAT Collaboration, F. Cafardo and R. Nemmen, “Fermi-LAT Observations of Sagittarius A*: Imaging Analysis,” Astrophys. J. 918 no. 1, (2021) 30, [arXiv:2107.00756 [astro-ph.HE]].
  * (56) CEPC Study Group Collaboration, M. Dong et al., “CEPC Conceptual Design Report: Volume 2 - Physics & Detector,” [arXiv:1811.10545 [hep-ex]].
  * (57) ILC Collaboration, “The International Linear Collider Technical Design Report - Volume 2: Physics,” [arXiv:1306.6352 [hep-ph]].
  * (58) FCC Collaboration, A. Abada et al., “FCC-ee: The Lepton Collider: Future Circular Collider Conceptual Design Report Volume 2,” Eur. Phys. J. ST 228 no. 2, (2019) 261–623.
  * (59) L. Linssen, A. Miyamoto, M. Stanitzki, and H. Weerts (editors), “Physics and Detectors at CLIC: CLIC Conceptual Design Report,” [arXiv:1202.5940 [physics.ins-det]].
  * (60) P. Ko and H. Yokoya, “Search for Higgs portal DM at the ILC,” JHEP 08 (2016) 109, [arXiv:1603.04737 [hep-ph]].
  * (61) Z. Heng, W. Wang, and H. Zhou, “Higgs-strahlung production of the lightest CP-even Higgs boson at ILC in natural NMSSM,” Chin. J. Phys. 55 (2017) 1723–1728.
  * (62) J. Liu, X.-P. Wang, and F. Yu, “A Tale of Two Portals: Testing Light, Hidden New Physics at Future $e^{+}e^{-}$ Colliders,” JHEP 06 (2017) 077, [arXiv:1704.00730 [hep-ph]].
  * (63) T. Kamon, P. Ko, and J. Li, “Characterizing Higgs portal dark matter models at the ILC,” Eur. Phys. J. C 77 no. 9, (2017) 652, [arXiv:1705.02149 [hep-ph]].
  * (64) D. Azevedo, M. Duch, B. Grzadkowski, D. Huang, M. Iglicki, and R. Santos, “Testing scalar versus vector dark matter,” Phys. Rev. D 99 no. 1, (2019) 015017, [arXiv:1808.01598 [hep-ph]].
  * (65) B. Grzadkowski, M. Iglicki, K. Mekala, and A. F. Zarnecki, “Dark-matter-spin effects at future $e^{+}e^{-}$ colliders,” JHEP 08 (2020) 052, [arXiv:2003.06719 [hep-ph]].
  * (66) G. Haghighat, M. Mohammadi Najafabadi, K. Sakurai, and W. Yin, “Probing a light dark sector at future lepton colliders via invisible decays of the SM-like and dark Higgs bosons,” Phys. Rev. D 107 no. 3, (2023) 035033, [arXiv:2209.07565 [hep-ph]].
  * (67) Y. Cheng, S.-F. Ge, J. Sheng, and T. T. Yanagida, “In preparation,”.