#  Probing LHAASO Galactic PeVatrons through gamma-ray and neutrino
correspondence

Prantik Sarmah1 Sovan Chakraborty,1 and Jagdish C. Joshi,2
1Department of Physics, Indian Institute of Technology Guwahati, Assam-781039,
India
2 Aryabhatta Research Institute of Observational Sciences (ARIES), Manora
Peak, Nainital 263001, India
prantik@iitg.ac.insovan<EMAIL_ADDRESS>
###### Abstract

Recently, Large High Altitude Air Shower Observatory (LHAASO) has detected
several Galactic point sources of ultra high energy (UHE; $E_{\gamma}>100$
TeV) gamma-rays. These gamma-rays are possibly created in leptonic or hadronic
interactions of cosmic rays (CRs) of PeV energies. In the hadronic channel
($p-p$ interaction), the gamma-rays are accompanied by neutrinos. The
detection of neutrinos is therefore crucial in understanding CR acceleration
in such objects. To estimate the neutrino flux, we adopt the two LHAASO
sources (J2226+6057, J1908+0621) found to be spatially associated with the
Supernova remnants (SNR G106.3+2.7, SNR G40.5-0.5). For these two sources, the
detected TeV-PeV gamma-ray spectra are found to be unusually hard (with
spectral index $\sim$ 1.8). We develop a model of gamma-ray and neutrino
emission based on the above two prototypes. The neutrino fluxes from these two
sources are found to be below the IceCube sensitivity, but are detectable in
upcoming IceCube-Gen2 and KM3NeT experiments. We further estimate the neutrino
fluxes from similar other 10 LHAASO PeVatron sources and most of them are
found to be detectable in IceCube-Gen2 and KM3NeT. Finally, we explore our
model parameters, in particular the spectral power law index and estimate the
future potential of the neutrino detectors to probe CR acceleration in such
Galactic sources.

###### keywords:

Cosmic rays (CRs) – Supernova remnant (SNR) – Gamma-rays – Neutrinos

††pubyear: 2015††pagerange:  Probing LHAASO Galactic PeVatrons through gamma-
ray and neutrino correspondence – Probing LHAASO Galactic PeVatrons through
gamma-ray and neutrino correspondence

##  1 Introduction

Recently, ground based Cherenkov telescope, Large High Altitude Air Shower
Observatory (LHAASO) has observed UHE (TeV-PeV) gamma-rays from 12 objects
located in our Galaxy (Abeysekara et al., 2020; Cao et al., 2021). One of
these sources is also detected by the Tibet AS+MD experiment and found to be
associated with the supernova remnant (SNR) G106.3+2.7 (Tibet AS$\gamma$
Collaboration et al., 2021). In past, other gamma-ray experiments like HAWC,
VERITAS and HESS have detected TeV gamma-rays from some of these sources
(Acciari et al., 2009; H. E. S. S. Collaboration et al., 2018; Abeysekara et
al., 2020). These gamma-rays are probably connected to the knee region of the
observed cosmic ray (CR) spectrum at around 3 PeV and are expected to be
produced by the Galactic sources (Gaisser et al., 2013). Supernova remnants
(SNRs) and pulsars are considered to be the most favourable Galactic objects
producing CRs around PeV energies (Bell, 2004; Reynolds, 2008; Gaggero et al.,
2018; Gabici et al., 2019; Cristofari, 2021). The LHAASO collaboration paper
(Cao et al., 2021), have listed the SNRs and pulsars that are spatially
associated with the observed gamma-rays and are consistent with the expected
potential sources of CRs (Montmerle, 1979; Gabici & Aharonian, 2007; Aharonian
et al., 2019). Individual supernova remnants (SNRs), i.e., IC 443 and W44, in
our Galaxy have also been detected by the Fermi-Satellite in the MeV-GeV band,
providing important hints of proton acceleration in SNRs (Ackermann et al.,
2013). Further, decade long observations of the Galactic centre region by the
HESS telescope have indicated that Sagittarius A∗ could accelerate protons
upto PeV energies (HESS Collaboration et al., 2016).

Another interesting source in our Galaxy is CRAB pulsar wind nebula (PWN), in
which the radiation has been observed from radio wavelengths to PeV gamma-
rays. LHAASO collaboration have also detected high energy photons from CRAB
PWN upto a maximum energy of $(1.12\pm 0.09)$ PeV (Lhaaso Collaboration et
al., 2021). They have interpreted the multi-wavelength radiation using a
synchrotron self-Compton model with electrons having maximum energy $\approx
2.15$ PeV, hence termed as leptonic PeVatron. However, the slightly larger
gamma-ray flux at the tail of the UHE spectra in CRAB PWN is considered to be
of hadronic origin (CR proton acceleration) (Liu & Wang, 2021; Peng et al.,
2022). Hence, CRAB PWN might be a source of both CR electrons and protons upto
PeV energies. Indeed, the observed CR proton data by Dark Matter Particle
Explorer (DAMPE) satellite indicates more than one class of Galactic CR
sources (An et al., 2019; Huang et al., 2020).

The gamma-rays detected by LHAASO can also have the possibility of both
leptonic or hadronic origin (Tibet AS$\gamma$ Collaboration et al., 2021; Cao
et al., 2021). In the leptonic channel, such gamma-rays may be produced by
inverse Compton of relativistic electrons. For the hadronic channel, the
interaction of CR protons with background protons in dense gas medium can
provide a dominant contribution. Gamma-rays and neutrinos are produced
simultaneously in this hadronic interaction (Gabici & Aharonian, 2007; Fujita
et al., 2009). Detection of the associated neutrinos together with the gamma-
rays will probe this hadronic origin (Gupta, 2013; Celli et al., 2020) and
will differentiate leptonic channel.

For the hadronic channel, interaction of CR protons accelerated in SNR with
molecular cloud is considered to be a prominent source of the UHE flux. SNRs
with age ($\rm t_{\rm age}\sim 100$ yr) are favoured as accelerators of CRs
upto PeV energies due to amplification of the magnetic field in the downstream
shock regions (Gaggero et al., 2018). The acceleration depends on the
confinement time of CR protons and gas density of the acceleration zone.
Longer confinement would provide the protons with sufficient time to
accelerate to higher energies. The association of UHE gamma-ray sources with
old SNRs ($\rm t_{\rm age}\geq 10$ kyr) might indicate a long time confinement
of CR protons, surrounding the SN shocks (Kar & Gupta, 2022; De Sarkar &
Gupta, 2022). The gas density determines magnetic field strength and also
responsible for CR energy losses due to interaction like $p-p$ collisions. For
example, for a confinement time of about 10 kyr, the required source gas
density is found to be about 10 $\rm cm^{-3}$ (see Sudoh & Beacom, 2022, for
details). These accelerated protons are expected to escape the source and
interact with neighbouring molecular clouds. Therefore, detection of neutrinos
produced in these interactions is crucial to probe the origin of these
PeVtrons.

These neutrinos might be detected by current (IceCube) and future (IceCube-
Gen2, KM3NeT) high energy neutrino detectors. In fact, it has been shown that
IceCube and KM3NeT are sensitive to objects like SNR G40.5-0.5 and Vela junior
(Mandelartz & Becker Tjus, 2015). IceCube has detected bunch of diffuse
neutrinos over a period of 7.5 years (Abbasi et al., 2021). Although, most of
these neutrinos are found to be of extra-galactic origin (Sarmah et al., 2022;
Chakraborty & Izaguirre, 2015; Petropoulou et al., 2017), the possibility of
few of these events originating in the galaxy can not be ruled out (Razzaque,
2013; Neronov et al., 2014). In addition, detection of gamma-rays by LHAASO
and TibeT AS+MD makes Galactic PeVatrons strong contenders of high energy
neutrino factories, detection of neutrinos can reveal these gamma-ray
associated PeVatrons.

In this work, we use the $p-p$ interaction model to estimate the resultant
gamma-ray and neutrino fluxes. We built our model in the context of the two
sources, SNR G106.3+2.7 and SNR G40.5-0.5. Both these sources have associated
molecular cloud in close vicinity. The gamma-rays mostly illuminated by the
molecular cloud motivates for a hadronic model with protons acceleration and
diffusing out of the SNR interacts with the protons in molecular cloud. We
find that neutrino flux from both these objects are below the IceCube
sensitivity. We also estimate the associated neutrino flux to observed gamma-
rays from the Crab nebula (Lhaaso Collaboration et al., 2021), considering
that the observed gamma-ray flux is completely hadronic. Even with this purely
hadronic origin, the Crab neutrino flux is below IceCube limit. We further
compute the associated neutrino fluxes from the remaining LHAASO sources based
on the observed gamma-ray flux. All these objects are also found to be below
IceCube sensitivity. However, many of these objects are found to be potential
neutrino sources for the future neutrino detectors like IceCube-Gen2 and
KM3NeT. We extended our study for the possible parameter space of our
molecular cloud model; mainly the spectral index and total energy.

This manuscript is organised as follows. In Sec. 2, we describe the model of
secondary (gamma-ray and neutrino) emission in the hadronic channel and fitted
the observed gamma-ray data of the two prototype sources (SNR G106.3+2.7 and
SNR G40.5-0.5). Sec. 3 is dedicated to discussion on the available parameter
space and detection prospects of Galactic PeVatrons with different gamma-ray
and neutrino telescopes. Finally, we conclude this paper in Sec. 4.

Figure 1: Hadronic gamma-rays (left), and muon neutrinos (right) from SNR
G106.3+2.7. The high gamma-ray flux is fitted to the Tibet AS+MD data for a
primary proton spectra with index 1.8 and the total energy of proton is
$8\times 10^{47}$ erg. The low energy part of the spectra is fitted to the
Fermi-LAT and VERITAS data. The absorption due to pair production on ISRF is
negligible that results in overlapping of the gamma-ray flux with (red-dashed)
and without absorption (blue). The distance to the source is taken to be 0.8
kpc. The correlated muon neutrino flux has been plotted (right) in addition to
the detectors’ sensitivities. The sensitivities of the neutrino detectors are
chosen for declination, $\delta=60^{\degree}$ for IceCube and
$\delta=30^{\degree}$ for IceCube-Gen2 which are closest available to the
declination of this SNR. IceCube is not sensitive to this neutrino flux but
IceCube-Gen2 has good sensitivity.

Figure 2: Hadronic gamma-rays (left), and muon neutrinos (right) from SNR
G40.5-0.5. The gamma-ray flux is fitted to the LHAASO J1908+0621 data for a
primary proton spectra with index 1.75 and the total energy of protons is
$9\times 10^{49}$ erg. The distance to the source is taken to be 8.5 kpc.
Gamma-rays above 10 TeV suffer tiny loss due to pair production on ISRF that
is shown by the red-dashed curve while the gamma-ray spectra without
absorption is shown by the blue curve. The corresponding neutrino flux has
good sensitivity to IceCubt-Gen2 but weak to IceCube. The sensitivities of the
neutrino detectors are shown for the declination, $\delta=0^{\degree}$ closest
available to the declination of this SNR.

##  2 Secondary Gamma-ray and Neutrino Fluxes

The interaction of CR protons with the molecular cloud leads to the production
of charged $\pi^{\pm}$ and neutral pions $\pi^{0}$, respectively. The decay of
neutral pions produces gamma-rays ($\gamma$) and charged pions produce
neutrino ($\nu$) fluxes. The fluxes of these gamma-rays and neutrinos depend
on the CR proton spectra ($J_{p}(E_{p})$), the molecular cloud density
($N_{H}$), and the production spectra of the secondary particles
($\gamma,~{}\nu$) for a given energy of proton ($E_{p}$) (Kelner et al.,
2006). The flux of secondary particles is expressed as

$\Phi_{i,pp}(E_{i})=cn_{H}\int_{\rm
E_{i}}^{\infty}\sigma_{pp}(E_{p})J_{p}(E_{p})F_{i}\left(\frac{E_{i}}{E_{p}},E_{p}\right)\frac{dE_{p}}{E_{p}},$
(1)

where, $i=\gamma~{}{\rm or}~{}\nu_{f}$ and $f$ stands for neutrino flavours
($e$, $\mu$). The CR protons are assumed to follow a power-law type
distribution with an exponential cutoff above $E_{0}$, i.e.,
$J_{p}(E_{p})=A_{p}E_{p}^{-\alpha}{\rm exp}(-E_{p}/E_{0})$ $\rm TeV^{-1}$. The
power-law index $\alpha$ is typically chosen to be 2 based on Fermi’s
diffusive shock acceleration mechanism (Gaisser et al., 2016). The
normalization constant $A_{p}$ is related to the total energy ($E_{p,\rm
total}$) of the CR protons by the relation
$\int_{m_{p}}^{\infty}E_{p}J_{p}(E_{p})\mathrm{d}E_{p}=E_{p,\rm total}$ and
$\sigma_{pp}(E_{p})$ is the energy dependent cross-section for the $p-p$
interaction (Kafexhiu et al., 2014). In particular, the upper limit of this
integration is the maximum proton energy, $E_{p,\rm max}$ which depends on the
magnetic field of the acceleration zone, confinement time and shock speed
(Gaisser et al., 2016; Tibet AS$\gamma$ Collaboration et al., 2021). In case
of SNR, this $E_{p,\rm max}$ is found to be a few PeV (Tibet AS$\gamma$
Collaboration et al., 2021). In our analysis, we obtain the proton
normalization by fitting the observed gamma-rays, and thus the exact value of
$E_{p,\rm max}$ is not necessary for our calculation. The cut-off energy
$E_{o}<E_{p,\rm max}$ makes the proton spectra to fall rapidly for energies
larger than $E_{o}$.

Now, these secondary gamma-rays and the neutrinos are produced from the same
$p-p$ interaction implying a connection between the final state fluxes of
these secondaries. The secondary fluxes are connected by the following
relation at $E_{\gamma}\approx 2E_{\nu}$ (Ahlers & Murase, 2014)

$E_{\gamma}\frac{\mathrm{d}N_{\gamma}(E_{\gamma})}{\mathrm{d}E_{\gamma}}\simeq
e^{-d/\lambda_{\gamma\gamma}}\frac{2}{3}\sum_{f}E_{\nu_{f}}\frac{\mathrm{d}N_{\nu_{f}}(E_{\nu})}{\mathrm{d}E_{\nu_{f}}},$
(2)

where, $\frac{\mathrm{d}N_{\gamma}(E_{\gamma})}{\mathrm{d}E_{\gamma}}$ and
$\frac{\mathrm{d}N_{\nu_{f}}(E_{\nu})}{\mathrm{d}E_{\nu_{f}}}$ are the
differential fluxes of gamma-rays and neutrinos of flavour $f$ respectively.
$\lambda_{\gamma\gamma}$ is the mean free path for the UHE gamma-rays in the
Galactic radiation fields for a source located at a distance $d$. Therefore,
one might predict the neutrino flux from SNR using this multi-messenger
connection, i.e., from gamma-ray observations.

The same secondary gamma-ray and neutrino flux at earth can also be estimated
directly from Eq. 1 and given by

$E_{i}^{2}\frac{\mathrm{d}N_{i}}{\mathrm{d}E_{i}}=\frac{E_{i}^{2}\Phi_{i,pp}(E_{i})}{4\pi
d^{2}}.$ (3)

Since we are interested in detection of muon tracks with high energy neutrino
detectors like IceCube, the muon flux at earth,
$E_{\nu_{\mu}}^{2}\frac{\mathrm{d}N_{\nu_{\mu}}}{\mathrm{d}E_{\nu_{\mu}}}$ is
computed by assuming the flavour ratio $1:1:1$. Note that neutrinos and anti-
neutrinos are treated as same here. Clearly, for sources with known properties
like $n_{H}$, $E_{0}$ etc., one may use Eq. 3, otherwise the Eq. 2 can be used
for sources with only gamma-ray observations.

Neutrinos being weakly interacting in nature can propagate to earth without
any losses. However, gamma-rays may undergo absorption during propagation.
Absorption might effect the gamma-rays produced in SNRs and therefore, it is
important to analyse this absorption process. This might also effect the cut-
off energy, $E_{0}$ of the CR protons i.e., large absorption would allow large
$E_{0}$ for the primary CR flux. The absorption is due to pair production on
low energy photon backgrounds like cosmic microwave background (CMB) and the
Galactic interstellar radiation field (ISRF). This is estimated by the factor
$e^{-\tau_{\gamma\gamma}}$, where
$\tau_{\gamma\gamma}=d/{\lambda_{\gamma\gamma}}$ is the attenuation of gamma-
rays on the low energy photons (Moskalenko et al., 2006). Larger the optical
depth, larger is the absorption. The optical depth depends on the density and
average energies of these low energy photons and the distance to the source
$d$. Clearly, gamma-ray fluxes from further sources experience larger
absorption. The CMB photons’ number density is about $440~{}\rm cm^{-3}$ and
has an average energy of the order of $10^{-4}$ eV. This results in threshold
energy for gamma-rays to interact with CMB photons at about 100 GeV
(Moskalenko et al., 2006). Since the CMB is uniform, the amount of absorption
on CMB effectively depends only on the distance, $d$. On the other hand, the
energies of Galactic ISRF photons are larger than that of CMB photons and the
distribution of ISRF photons is highly non-uniform. Therefore, the threshold
energy for ISRF is lower than 100 GeV and the absorption depends on the
location of the source (Moskalenko et al., 2006). For example, for a source
located in a denser ISRF region will have larger absorption than that of a
source located in a thinner ISRF.

In the following, we use the method of direct calculation, i.e., Eqs. 1 and 3
for the secondary fluxes from the Galactic SNR sources, G106.3+2.7 and
G40.5-0.5, as requisite properties are known from observations.

### SNR G106.3+2.7

SNR G106.3+2.7 is one of the probable PeVatron source in which CR protons are
accelerated due to the high shock speed (Ge et al., 2021). The distance of
this source is 0.8 kpc and its shape is not perfectly spherical but a comet
shaped SNR of length 14 pc and width of 6 pc (Kothes et al., 2001). Recently,
The Tibet AS$\gamma$ Collaboration reported the discovery of gamma-rays above
10 TeV which extends upto 100 TeV gamma-rays (Tibet AS$\gamma$ Collaboration
et al., 2021). Further, the LHAASO collaboration also reported the discovery
of UHE gamma-rays from this spatial region from the source LHAASO J2226+6057
(Cao et al., 2021). These UHE gamma-rays are found to be spatially correlated
with a molecular cloud and hence CR proton interactions are considered to be
strong contenders (Tibet AS$\gamma$ Collaboration et al., 2021). This is also
supported by the hadronic interaction models developed for this source (Bao &
Chen, 2021; Yang et al., 2022). The probable density of the molecular cloud
region is approximately $10~{}{\rm cm^{-3}}$ and is consistent with the
Bremsstrahlung radiation of electrons (Fujita et al., 2021).

Figure 1 shows the model fitting of Tibet AS+MD gamma-ray data (left) and the
corresponding neutrino emission in the model (right). The multi-messenger
nature of the emission is crucial in fitting the gamma-ray spectra. The
spectra data at the lower energies ($10^{-3}-1$ TeV) are from VERITAS (Acciari
et al., 2009), at intermediate energies ($1-10$ TeV) are from Fermi-LAT (Xin
et al., 2019), and at higher energies ($\sim 10$ TeV) are taken from Tibet
AS$\gamma$ (Tibet AS$\gamma$ Collaboration et al., 2021). The best fitting
parameters $E_{p,\rm total}$, $E_{0}$ and $\alpha$ are such that both the
normalisation and the slope of the spectra are well explained (see Table 1).
We note the smaller $E_{p,\rm total}$ for this SNR and this may be due to low
transfer of shock kinetic energy to accelerated protons (Tibet AS$\gamma$
Collaboration et al., 2021). Moreover, the power law index, $\alpha$ in this
object is found to be less than 2 from gamma-ray spectra fitting. This is in
contrast with the diffusive shock acceleration theory that predicts
$\alpha=2$. Extreme pressure at shock front might cause such efficient
particle acceleration (see e.g., Malkov, 1999a). These fit parameters are also
consistent with (Tibet AS$\gamma$ Collaboration et al., 2021). The source
being located only at 0.8 kpc, the gamma-ray flux has negligible absorption
during propagation to earth (Tibet AS$\gamma$ Collaboration et al., 2021;
Moskalenko et al., 2006).

Based on these parameter choices, we estimate the gamma-ray flux from this
object. The resulting gamma-ray flux gives a good fitting to the gamma-ray
data at both low energy and high energy. As mentioned, absorption during
propagation is negligible which resulted in the overlap of the blue and red
dashed curves. The associated neutrino flux is estimated using Eq. 3 and is
found to be consistent with Eq. 2. This neutrino flux is found to be much
below the detection sensitivity of IceCube. However, IceCube-Gen2 will be able
to detect these neutrinos provided the source is located closer to earth. Note
that larger $\alpha$ would change the sensitivity.

### SNR G40.5-0.5

SNR G40.5-0.5 is spatially overlapped with the UHE gamma-ray source LHAASO
J1908+0621 (Cao et al., 2021). In earlier observations, this SNR was found be
overlapped with extended TeV source MGRO J1908+06, detected by Milagro
Galactic plane survey (Abdo et al., 2007). This source was also detected by
the H.E.S.S. telescope in gamma-rays (Aharonian et al., 2009). UHE gamma-rays
were also discovered from MGRO J1908+06 source by the HAWC detector
(Abeysekara et al., 2020) that makes SNR G40.5-0.5 a plausible CR particle
accelerator. The source distance is uncertain and CO Molecular-line emission
infers a source distance of approximately 3.4 kpc (Yang et al., 2006) but at a
larger distance of 5.5 to 8.5 kpc using $\Sigma-$D relation (Downes et al.,
1980). The mean gas density in the associated molecular cloud is $\sim
45~{}{\rm cm^{-3}}$ taking the source distance 8 kpc (Li et al., 2021) and the
source linear size at this distance would be approximately 60 pc.

Model fit of the gamma-ray data of this object is shown in the left panel Fig.
2 for the parameters listed in Table 1. The low energy ($10^{-3}-10^{-1}$ TeV)
data are from Fermi-LAT (Abdollahi et al., 2020; Cao et al., 2021) and at high
energies ($>10^{-1}$ TeV) data are taken from HESS (H. E. S. S. Collaboration
et al., 2018), HAWC (Abeysekara et al., 2020), and LHAASO (Cao et al., 2021).
The parameters listed in Table 1 can well explain these multi-messenger gamma-
ray data. Note that $E_{p,\rm total}$ in this SNR is found to be larger than
that of SNR G106.3+2.7. However, the power-law index $\alpha$ is found to be
similar to that of G106.3+2.7. Since this source is farther away, gamma-rays
suffer larger absorption due to pair production on ISRF compared to
G106.3+2.7. However, this absorption is still tiny and negligible; and shown
by the red dashed curve. The corresponding neutrino flux is estimated using
Eq. 3 and shown in the right panel. This neutrino flux estimation is found to
be consistent with Eq. 2. To check the IceCube neutrino sensitivity, we look
into the closest available declination of IceCube, i.e., $\delta=0^{\degree}$
to the source declination. The neutrino flux has weak sensitivity to IceCube
but has reasonably good sensitivity to IceCube-Gen2. Hence, this object sets
an outstanding goal to neutrino detection in IceCube-Gen2 as the neutrino flux
is much above the detection sensitivity.

In the following, we investigate the parameter space of Galactic PeVatrons
allowed by gamma-ray observations (Tibet AS$\gamma$ Collaboration et al.,
2021; Cao et al., 2021) and discuss the future detection prospects with gamma-
ray and neutrino detectors.

Parameters | SNR G106.3+2.7 | SNR G40.5-0.5
---|---|---
$E_{p,\rm total}$ (erg) | $8\times 10^{47}$ | $9\times 10^{49}$
$E_{{}_{0}}$ (TeV) | $6\times 10^{2}$ | $4\times 10^{2}$
$\alpha$ | 1.8 | 1.75
Table 1: Model parameters chosen for the two prototype SNRs. These parameters
give the best fit to the observed gamma-ray data of these objects. SNR
G106.3+2.7 is located at 0.8 kpc and G40.5-0.5 is located at 8.5 kpc. The
target gas density, $n_{H}$ for SNR G106.3+2.7 is $10~{}\rm cm^{-3}$ and for
G40.5-0.5 is $45~{}\rm cm^{-3}$.

Figure 3: Left: Detection prospects of gamma-rays from Galactic SNR at 10 kpc
with different gamma-ray telescopes. The curves in different colours shows the
sensitivities of Fermi-LAT (green) (Fermi-LAT, ), LHAASO (magenta) (Vernetto &
for the LHAASO Collaboration, 2016), Tibet AS+MD (brown) (Tibet, ) and CTA
(black) (CTA, ). The orange band shows gamma-ray flux for the parameter space
based on the observations of SNR G106.3+2.7 and G40.5-0.5. This shows that the
high energy gamma-ray detectors will probe large fraction of this parameter
space for SNR at 10 kpc. Right: Corresponding detection prospects at high
energy neutrino detectors. The sensitivities of IceCube, IceCube-Gen2 and
KM3NeT are shown by the blue dot-dashed, magenta and green lines respectively
(IceCube Collaboration et al., 2021; Aartsen et al., 2021; Aiello et al.,
2019). These neutrino detectors will also probe some part of the parameter
space.

Figure 4:  Left: Detection horizon of SNRs for different gamma-ray telescopes.
The orange band show integrated flux above 1 TeV as a function of distance
(kpc) for the same parameter as in Fig. 3.. The horizontal lines show the
detectors sensitivities (Vernetto & for the LHAASO Collaboration, 2016; Fermi-
LAT, ; CTA, ; Tibet, ). All these gamma-ray telescopes are capable of probing
a large part of the parameter space . The upcoming CTA will probe sources as
far as 50 kpc. Right: Corresponding detection horizon of different neutrino
detectors. The orange band shows the energy integrated ($10^{1}-10^{2}$ TeV)
$\nu_{\mu}$ flux as a function of distance (kpc) for the same parameter as in
Fig. 3. The horizontal lines show the energy integrated sensitivities
(corresponds to $\delta=0^{\degree}$ in the energy range $10^{1}-10^{2}$ TeV)
of IceCube (blue), IceCube-Gen2 (magenta) and KM3NeT (green) (IceCube
Collaboration et al., 2021; Aartsen et al., 2021; Aiello et al., 2019). This
shows that the maximum reach of IceCube is about 10 kpc whereas for IceCube-
Gen2 and KM3NeT, the maximum reach is about 50kpc. The stars of different
colours in both panels show the energy integrated $\gamma$ and $\nu_{\mu}$
flux for the sources listed in Table 2.

##  3 Detection prospects of Galactic PeVatrons

In the preceding sections, we have computed the expected muon neutrino flux in
detail for the two prototype sources. The neutrino flux estimates of such
objects show that the flux might be detected by the current and upcoming
generation of IceCube experiment. This motivates one to look further into the
detection prospects of these Galactic objects. Based on this, we make estimate
of future detection of such sources motivated by the two detected prototype
SNRs. We also estimate the expected neutrino events in IceCube from the LHAASO
catalog sources (see Table 2 and (Cao et al., 2021)) to check the consistency
of our chosen parameter space.

In order to investigate the detection possibility, a range of CR proton
spectra parameters is considered, motivated by the observations (Cao et al.,
2021; Tibet AS$\gamma$ Collaboration et al., 2021). This parameter space is
obtained by varying the total energy $E_{p,\rm total}$ ($10^{48}-10^{51}$ erg)
and the power law index $\alpha$ ($1.7-2.0$) based on the gamma-ray
observations of SNR G106.3+2.7 and G40.5-0.5. While varying the parameters, we
choose to keep the cut-off energy, $E_{0}$ fixed as both objects show similar
values. The power law index, $\alpha$ is also found to be similar in these two
objects, i.e., about 1.8. However, this index is in contrast with the
diffusive shock acceleration theory which predicts a power law index
$\alpha=2.0$ (Gaisser et al., 2016). Such smaller value of $\alpha$ is
possibly due to very efficient proton acceleration in these objects (Malkov,
1999a; Blasi, 2002). Indeed, observation of neutrinos from such sources by the
neutrino detectors in question will be crucial to probe the softer spectra
i.e., large $\alpha$. Therefore, the variation of $\alpha$ is considered in
the range $(1.7-2.0)$.

For the parameter space discussed above, gamma-ray and muon neutrino fluxes
are estimated for an SNR at 10kpc. The left panel in Fig. 3 shows the gamma-
ray flux and the right panel shows the corresponding muon neutrino flux. The
orange band in both panels corresponds to the uncertainty in the parameters
$E_{p,\rm total}$ and $\alpha$. In the gamma-ray plot, the sensitivities of
different gamma-ray telescopes operating at different energy ranges have been
shown. At low energy, Fermi-LAT is capable of probing a small part of the SNR
parameter space. Interestingly, the high energy gamma-ray telescopes like
LHAASO, Tibet AS+MD and CTA (upcoming) will be able to probe a much larger
fraction of the SNR parameter space.

In addition to the correlated neutrino flux, high energy neutrino detectors’
sensitivities have also been plotted in right panel of Fig. 3. The blue dot
dashed line corresponds to the sensitivity of IceCube (IceCube Collaboration
et al., 2021), the magenta line shows the sensitivity of IceCube-Gen2 (Aartsen
et al., 2021) and the green-dashed line is the sensitivity of KM3NeT (Aiello
et al., 2019). This shows that IceCube is hardly sensitive to such neutrinos
from sources at 10 kpc. However, IceCube-Gen2 and KM3NeT might be able to
detect these neutrinos from SNRs provided the neutrino fluxes produced are
$>(10^{-13})\rm~{}TeVcm^{-2}s^{-1}$ above 10 TeV. Nevertheless, these future
detectors together with the gamma-ray telescopes can put stringent constraints
on a significant region of the SNR parameter space. For this, the UHE gamma-
rays detected by LHAASO (Cao et al., 2021) might assist to scan this parameter
space and help us to understand the future detection possibilities.

Sl No. | LHAASO sources | IceCube events, $\mathcal{N}_{\nu_{\mu}}$ (5 yr)
---|---|---
1 | J0534+2202 | 2.74
2 | J1825-1326 | 0.28
3 | J1839-0545 | 0.08
4 | J1843-0338 | 4.68
5 | J1849-0003 | 2.82
6 | J1908+0621 | 4.55
7 | J1929+1745 | 0.85
8 | J1956+2845 | 1.33
9 | J2018+3651 | 2.33
10 | J2108+5157 | 1.05
11 | J2032+4102 | 1.07
12 | J2226+6057 | 2.57
Table 2: Expected number $\nu_{\mu}$ events at IceCube for the LHAASO
catalogue. The sources J2226+6057 and J1908+0621 are associated with SNR
G106.3+2.7 and G40.5-0.5 respectively. The conventional atmospheric
$\nu_{\mu}$ flux produces a background of $\mathcal{O}(10^{4})$ events on top
of these events.

The PeVatron sources detected in gamma-rays by LHAASO (Cao et al., 2021) are
listed in Table 2. The estimation of neutrinos from these sources is not
straightforward due to the following reasons. All of these objects might not
be pure hadronic PeVatrons and therefore might not produce much neutrinos
(Sudoh & Beacom, 2022). In fact, the source J0534+2202 is found to be
associated with the Crab nebula and most of these gamma-rays are expected to
be of leptonic origin except for the tail of the spectra which might be
hadronic (Lhaaso Collaboration et al., 2021). Therefore, we use the gamma-ray-
neutrino correlation given by Eq. 2 to estimate the maximum possible neutrinos
from the Crab considering these gamma-rays to be completely hadronic. For the
remaining LHAASO sources, association of any molecular cloud is not well known
(Cao et al., 2021) and model of CR interacting with molecular cloud might not
be appropriate. In addition, only the fluxes at 100 TeV are available for
these sources (Lhaaso Collaboration et al., 2021). Thus, we assume a proton
power law index $\alpha=1.8$ and normalize the gamma-ray flux at 100 TeV to
estimate the neutrino fluxes from these objects (Cao et al., 2021).

To analyze the detection possibility of these sources with IceCube, we have
computed the number of $\nu_{\mu}$ tracks at IceCube by the following formula

$\mathcal{N}_{\nu_{\mu}}=T\int_{0.1~{}TeV}^{100~{}TeV}\mathrm{d}E_{\nu_{\mu}}\frac{\mathrm{d}N_{\nu_{\mu}}}{\mathrm{d}E_{\nu_{\mu}}}A_{eff}(\delta)\,$
(4)

where, $T$ is the observation time and $A_{eff}(\delta)$ is the declination
dependent effective area (Aartsen et al., 2017). The declination of the source
plays an important role in the detection of neutrinos. Therefore, the
declination of each source in Table 2 is taken into account for the
calculation of IceCube events. The number of events for 5 years are shown in
Table 2. The events are found to be so tiny that they might be buried deep in
the large conventional atmospheric background ($\mathcal{O}(10^{4})$). This
explains the non-observation of neutrinos till date from such SNRs (Abbasi et
al., 2022).

It is clear from the above analysis that the source distance is also an
important parameter for the detection of these PeVatron sources by gamma-ray
and neutrino telescopes. Fig. 4 shows the energy integrated gamma-ray (left)
and muon neutrino (right) flux as a function of source distance (in kpc). The
orange band in both panels corresponds to the uncertainty in the SNR
parameters, $E_{p,\rm total}$ and $\alpha$ as discussed above. The gamma-ray
flux is integrated above 1 TeV following the best sensitive energy ranges of
the different detectors. The horizontal lines in the left plot are the
integral sensitivities for CTA, LHAASO and Tibet AS+MD. All these gamma-ray
telescopes are well capable of detecting SNRs in the galaxy. In future, CTA
will probe vast region of the SNR parameter space.

Similarly, detection horizons of different high energy neutrino telescopes
like IceCube, IceCube-Gen2 and KM3NeT are shown in the right panel of Fig. 4.
The orange band shows the integrated $\nu_{\mu}$ flux in the range
($10-10^{2}$) TeV for the parameter space discussed above. Energy integrated
sensitivities in the range $10^{1}-10^{2}$ TeV of IceCube, IceCube-Gen2 and
KM3NeT for declination angle $\delta=0^{\degree}$ have also been plotted. This
plot shows that the maximum reach of IceCube is about 10 kpc whereas for
IceCube-Gen2 and KM3NeT the maximum reach is about 50kpc. However, for the
lower limit of the neutrino flux, the detection horizon for IceCube falls
below 200 pc and for IceCube-Gen2, KM3NeT, it falls below 1 kpc.

The gamma-ray and neutrino fluxes produced by the LHAASO sources (Cao et al.,
2021) listed in Table. 2 are also shown in Fig. 4 by the star symbols in
different colours. All these sources are above the sensitivities of LHAASO and
CTA, while they are below the sensitivity of IceCube. Interestingly, the
upcoming neutrino detectors IceCube-Gen2 and KM3NeT will be able to detect
neutrinos from such SNRs. As pointed out earlier, all of these sources might
not be pure hadronic PeVatrons (Sudoh & Beacom, 2022). Therefore, the neutrino
fluxes shown here can only be treated as upper limits. As mentioned earlier,
the Crab gamma-ray fluxes are expected to be mostly leptonic (Lhaaso
Collaboration et al., 2021). Thus, neutrino emission from Crab might be faint
and not detectable. Even for the purely hadronic model, the neutrino flux from
Crab is found to be below IceCube sensitivity. Indeed, IceCube-Gen2 and KM3NeT
will be able to constrain the hadronic nature of these sources.

Note that the detection of neutrinos strongly depends on the source
declination, $\delta$. The sensitivities shown here for the neutrino detectors
are for $\delta=0^{\degree}$, i.e., the largest sensitivities. For declination
angle other than this, the sensitivity would be lower and some of the sources
might fall below the detectors’ sensitivities. Nevertheless, the upcoming
neutrino detectors might still be able to detect some PeVatrons and unravel
the nature of these sources, i.e, hadronic or leptonic.

##  4 Discussions and Conclusions

The 12 sources detected by LHAASO in TeV gamma-rays provide evidence of CR
acceleration to very high energies. If these gamma-rays are created by
hadronic interaction of CR protons with background protons, they should be
accompanied by neutrinos. However, no neutrinos are expected for leptonic
origin of these gamma-rays. Therefore, detection of neutrinos will play a
crucial role to identify the nature of these sources. Based on the gamma ray
data, we have estimated the plausible neutrino fluxes from these sources and
test their detectability in IceCube, IceCube-Gen2 and KM3NeT.

These objects detected by LHAASO are found to be spatially associated with
SNRs and pulsars. CRs accelerated in SNRs and pulsar interacting with any
nearby molecular cloud can be the origin of these gamma-rays. In fact, for two
of these sources, gamma-rays are found to be originated in the molecular
clouds of SNR G106.3+2.7 and G40.5-0.5. However, for some of these sources no
molecular clouds have been identified yet. On the other hand, one of these
objects is found to be connected to the Crab nebula where the emission is
expected to be mostly leptonic. Therefore, detection of the associated
neutrino signal is important to probe the nature of these sources.

Using the gamma-ray data of SNR G106.3+2.7 and SNR G40.5-0.5, we have
developed a $p-p$ interaction model and estimated the neutrino fluxes from
these objects. The age of SNR G106.3+2.7 is $\sim 10$ kyr, while for SNR
G40.5-0.5, it is in between 10-20 kyr (see e.g., Cao et al., 2021). For these
two sources, model parameters, i.e., CR proton energy $E_{\rm p,\rm total}$,
spectral index $\alpha$, the cutoff energy $E_{0}$ are listed in Table 1. We
found $\alpha$ and $E_{0}$ to be similar for both SNRs (1.8 and $\sim 500$
TeV) but their $E_{\rm p,\rm total}$ is very different. The asymmetry in the
source geometry and cloud morphology can affect the escape of the CR protons
and may cause this difference in the total energy in CR protons (Fujita et
al., 2021). Our chosen spectral indices are harder compared to the standard
diffusive shock acceleration model, which predicts $\alpha=2$ (Drury, 1983).
However, if the shock compression ratio is very high then the CR spectral
index can be harder; i.e. $\alpha=1.5$ (Malkov, 1999b) as also predicted using
nonlinear diffusive shock acceleration mechanism (Blasi, 2002). The neutrino
flux from SNR G106.3+2.7 is found to be below IceCube sensitivity and SNR
G40.5-0.5 has poor sensitivity to IceCube. This also explains the non-
observation of neutrinos from these two sources. However, both these sources
appeared to be sensitive to IceCube-Gen2. Hence, IceCube-Gen2 will be able to
probe the hadronic channel in these objects.

Based on the analysis of these two sources, we have derived the probable
parameter space of this SNR + molecular cloud model. The parameter space has
been constructed considering uncertainties in the model parameters, i.e.,
spectral index, $\alpha$ and total energy, $E_{p,\rm total}$. The upper limit
of $\alpha=2.0$ is considered based on the standard CR acceleration mechanism.
The wide variation in $E_{p,\rm total}$ is taken based on the observation the
two SNR G106.3+2.7 and G40.5-0.5. Future neutrino telescopes (IceCube-Gen2 and
KM3NeT) will be able to probe a large portion of this parameter space,
provided the events are located within 50 kpc. However, the best possibility
of detection indeed lies within a region of 10 kpc which is inferred by the
LHAASO sources.

We have also estimated the neutrino flux from Crab nebula based on the gamma-
ray spectra observed by LHAASO. Most of these gamma-rays are expected to be of
leptonic origin. Thus, the neutrino signal is expected to be weak. The
neutrino flux from Crab is found to be below IceCube sensitivity, even for the
purely hadronic model. Neutrino fluxes for the remaining LHAASO sources are
estimated based on the gamma-ray data and assuming a power law index 1.8. All
these sources are found to be below IceCube sensitivity but many of them might
be detectable in IceCube-Gen2 and KM3NeT.

Further, primary and secondary electrons accelerated at the shock regions, can
also contribute towards the multi-wavelength radiation. Their synchrotron
radiation depends on the magnetic field and hence this can be tuned to get the
radio to X-ray spectrum consistent with the observations. The inverse-Compton
radiation can be another channel for the gamma-ray spectrum in the UHE regime.
In this work we are mainly interested in estimating the maximum neutrino flux
from SNRs, hence we have not discussed the radiation due to electrons.
However, contamination from the leptonic channel of the concentric PWN cannot
be discarded in these sources (Torres et al., 2014; Yu et al., 2022; Joshi et
al., 2022)

To conclude, we have analysed the detection prospects of the LHAASO sources
and also similar possible sources in the Galaxy. IceCube might be able to
detect closer objects within 1 kpc but for extremely neutrino bright sources
the detection horizon extends up to 10 kpc. The future detectors KM3NeT and
IceCube-Gen2 can have detection sensitivity up to 50 kpc. Such detection will
probe the nature of Galactic PeVatrons.

## Acknowledgements

S.C. acknowledges the support of the Max Planck India Mobility Grant from the
Max Planck Society, supporting the visit and stay at MPP during the project.
S.C has also received funding from DST/SERB projects CRG/2021/002961 and
MTR/2021/000540.

## References

  * Aartsen et al. (2017) Aartsen M. G., et al., 2017, The Astrophysical Journal, 835, 151
  * Aartsen et al. (2021) Aartsen M. G., et al., 2021, J. Phys. G, 48, 060501
  * Abbasi et al. (2021) Abbasi R., et al., 2021, Phys. Rev. D, 104, 022002
  * Abbasi et al. (2022) Abbasi R., et al., 2022, arXiv e-prints, p. arXiv:2211.14184
  * Abdo et al. (2007) Abdo A. A., et al., 2007, apjl, 664, L91
  * Abdollahi et al. (2020) Abdollahi S., et al., 2020, ApJS, 247, 33
  * Abeysekara et al. (2020) Abeysekara A. U., et al., 2020, Phys. Rev. Lett., 124, 021102
  * Acciari et al. (2009) Acciari V. A., et al., 2009, ApJ, 703, L6
  * Ackermann et al. (2013) Ackermann M., et al., 2013, Science, 339, 807
  * Aharonian et al. (2009) Aharonian F., et al., 2009, aap, 499, 723
  * Aharonian et al. (2019) Aharonian F., Yang R., de Oña Wilhelmi E., 2019, Nature Astronomy, 3, 561
  * Ahlers & Murase (2014) Ahlers M., Murase K., 2014, Phys. Rev. D, 90, 023010
  * Aiello et al. (2019) Aiello S., et al., 2019, Astropart. Phys., 111, 100
  * An et al. (2019) An Q., et al., 2019, Science Advances, 5, eaax3793
  * Bao & Chen (2021) Bao Y., Chen Y., 2021, apj, 919, 32
  * Bell (2004) Bell A. R., 2004, MNRAS, 353, 550
  * Blasi (2002) Blasi P., 2002, Astroparticle Physics, 16, 429
  * (18) CTA, Cherenkov Telescope Array, Available at https://www.cta-observatory.org/science/ctao-performance/#1472563157332-1ef9e83d-426c
  * Cao et al. (2021) Cao Z., et al., 2021, Nature, 594, 33
  * Celli et al. (2020) Celli S., Aharonian F., Gabici S., 2020, apj, 903, 61
  * Chakraborty & Izaguirre (2015) Chakraborty S., Izaguirre I., 2015, Phys. Lett. B, 745, 35
  * Cristofari (2021) Cristofari P., 2021, Universe, 7, 324
  * De Sarkar & Gupta (2022) De Sarkar A., Gupta N., 2022, ApJ, 934, 118
  * Downes et al. (1980) Downes A. J. B., Pauls T., Salter C. J., 1980, aap, 92, 47
  * Drury (1983) Drury L. O., 1983, Reports on Progress in Physics, 46, 973
  * (26) Fermi-LAT, Official website, Available at https://fermi.gsfc.nasa.gov
  * Fujita et al. (2009) Fujita Y., Ohira Y., Tanaka S. J., Takahara F., 2009, ApJ, 707, L179
  * Fujita et al. (2021) Fujita Y., Bamba A., Nobukawa K. K., Matsumoto H., 2021, apj, 912, 133
  * Gabici & Aharonian (2007) Gabici S., Aharonian F. A., 2007, ApJ, 665, L131
  * Gabici et al. (2019) Gabici S., Evoli C., Gaggero D., Lipari P., Mertsch P., Orlando E., Strong A., Vittino A., 2019, International Journal of Modern Physics D, 28, 1930022
  * Gaggero et al. (2018) Gaggero D., Zandanel F., Cristofari P., Gabici S., 2018, mnras, 475, 5237
  * Gaisser et al. (2013) Gaisser T. K., Stanev T., Tilav S., 2013, Frontiers of Physics, 8, 748
  * Gaisser et al. (2016) Gaisser T. K., Engel R., Resconi E., 2016, Cosmic Rays and Particle Physics, 2 edn. Cambridge University Press, doi:10.1017/CBO9781139192194
  * Ge et al. (2021) Ge C., Liu R.-Y., Niu S., Chen Y., Wang X.-Y., 2021, The Innovation, 2, 100118
  * Gupta (2013) Gupta N., 2013, Astroparticle Physics, 48, 75
  * H. E. S. S. Collaboration et al. (2018) H. E. S. S. Collaboration et al., 2018, A&A, 612, A1
  * HESS Collaboration et al. (2016) HESS Collaboration et al., 2016, Nature, 531, 476
  * Huang et al. (2020) Huang Z.-Q., Liu R.-Y., Joshi J. C., Wang X.-Y., 2020, apj, 895, 53
  * IceCube Collaboration et al. (2021) IceCube Collaboration et al., 2021, arXiv e-prints, p. arXiv:2101.09836
  * Joshi et al. (2022) Joshi J. C., Tanaka S. J., Miranda L. S., Razzaque S., 2022, arXiv e-prints, p. arXiv:2205.00521
  * Kafexhiu et al. (2014) Kafexhiu E., Aharonian F., Taylor A. M., Vila G. S., 2014, Phys. Rev. D, 90, 123014
  * Kar & Gupta (2022) Kar A., Gupta N., 2022, ApJ, 926, 110
  * Kelner et al. (2006) Kelner S. R., Aharonian F. A., Bugayov V. V., 2006, Phys. Rev. D, 74, 034018
  * Kothes et al. (2001) Kothes R., Uyaniker B., Pineault S., 2001, apj, 560, 236
  * Lhaaso Collaboration et al. (2021) Lhaaso Collaboration et al., 2021, Science, 373, 425
  * Li et al. (2021) Li J., et al., 2021, apjl, 913, L33
  * Liu & Wang (2021) Liu R.-Y., Wang X.-Y., 2021, ApJ, 922, 221
  * Malkov (1999a) Malkov M. A., 1999a, ApJ, 511, L53
  * Malkov (1999b) Malkov M. A., 1999b, apjl, 511, L53
  * Mandelartz & Becker Tjus (2015) Mandelartz M., Becker Tjus J., 2015, Astroparticle Physics, 65, 80
  * Montmerle (1979) Montmerle T., 1979, apj, 231, 95
  * Moskalenko et al. (2006) Moskalenko I. V., Porter T. A., Strong A. W., 2006, Astrophys. J. Lett., 640, L155
  * Neronov et al. (2014) Neronov A., Semikoz D., Tchernin C., 2014, Phys. Rev. D, 89, 103002
  * Peng et al. (2022) Peng Q.-Y., Bao B.-W., Lu F.-W., Zhang L., 2022, ApJ, 926, 7
  * Petropoulou et al. (2017) Petropoulou M., Coenders S., Vasilopoulos G., Kamble A., Sironi L., 2017, Mon. Not. Roy. Astron. Soc., 470, 1881
  * Razzaque (2013) Razzaque S., 2013, Phys. Rev. D, 88, 081302
  * Reynolds (2008) Reynolds S. P., 2008, araa, 46, 89
  * Sarmah et al. (2022) Sarmah P., Chakraborty S., Tamborra I., Auchettl K., 2022, JCAP, 08, 011
  * Sudoh & Beacom (2022) Sudoh T., Beacom J. F., 2022, arXiv e-prints, p. arXiv:2209.03970
  * (60) Tibet A., , The Tibet ASg Experiment, Available at icrr.u-tokyo.ac.jp/prwps/extreview/2013/pdf/16_04_Tibet_AS-Gamma.pdf
  * Tibet AS$\gamma$ Collaboration et al. (2021) Tibet AS$\gamma$ Collaboration et al., 2021, Nature Astronomy, 5, 460
  * Torres et al. (2014) Torres D. F., Cillis A., Martín J., de Oña Wilhelmi E., 2014, Journal of High Energy Astrophysics, 1, 31
  * Vernetto & for the LHAASO Collaboration (2016) Vernetto S., for the LHAASO Collaboration 2016, Journal of Physics: Conference Series, 718, 052043
  * Xin et al. (2019) Xin Y., Zeng H., Liu S., Fan Y., Wei D., 2019, ApJ, 885, 162
  * Yang et al. (2006) Yang J., Zhang J.-L., Cai Z.-Y., Lu D.-R., Tan Y.-H., 2006, cjaa, 6, 210
  * Yang et al. (2022) Yang C., Zeng H., Bao B., Zhang L., 2022, aap, 658, A60
  * Yu et al. (2022) Yu H., Wu K., Wen L., Fang J., 2022, na, 90, 101669


# Grid-Forming Control of Power Converters: Equivalence Proof through
Simplified Models

Aidar Zhetessov  Dept. of Electrical and Computer Engineering
1415 Engineering Drive
Madison, WI 53706

###### Abstract

This work establishes the equivalence of selected grid-forming control
algorithms within the context of simplified theoretical models. Considered
algorithms are droop control, Virtual Synchronous Machine (VSM) and matching
control. It is shown that nodal and network dynamics under those regulators
boil down to the same equations near the selected (trivial) nominal operating
point. Finally, some practical insights on each regulator dynamics and an
outlook are provided.

###### Index Terms:

grid-forming control, droop, VSM, matching, equivalence

##  I Introduction

Sustainable development pushes the technology towards grid integration of
renewable energy sources (RES) [1]. With more RES penetrating the grid, some
form of grid formation/support is required from the RES interfacing converters
for a stable grid operation [2]. As a response to that necessity, several
converter control algorithms have been developed with a general term “grid-
forming” control serving as an umbrella definition for those [3], [4].
Naturally arising questions here are how do those grid-forming techniques
compare and which one to use? Numerous detailed studies have addressed those
(e.g. [5]).

The small project at hand aims at contributing to the aforementioned works by
theoretically establishing the equivalence between the grid-forming control
techniques and gaining some practical insights on controller selection/tuning.
Considered techniques are droop, Virtual Synchronous Machine (VSM) and
matching control. To simplify the analysis, network AC power flows are
linearized around the nominal (trivial) network operating point. Moreover,
inner voltage-current regulators of the converters are assumed to be ideal.

The work is organized as follows: Section II introduces nodal dynamics with
grid-forming regulators and establishes equivalence. Section III analyzes
network dynamics and presents some insights on grid-forming regulator
selection/tuning from a network perspective. Finally, Section IV summarizes
the findings and shapes an outlook.

##  II Nodal Equivalence of Grid-Forming Controls

This section presents the nodal dynamics of the power converters under three
different grid-forming controllers: VSM, droop and matching. VSM is referred
to as reference regulator. Unless otherwise stated, all node dynamics are in
dq-frame (rotating at nominal grid frequency $\omega_{0}$), not linearized
(full) and virtual (state integration happens numerically).

###  II-A Virtual Synchronous Machine

VSM nodal dynamics are summarized in (1) for an arbitrary node $k$ of the
network, represented as a graph:

$\begin{split}\frac{d}{dt}\theta_{k}&=\omega_{k}\\\
M_{k}\frac{d}{dt}\omega_{k}&=-D_{k}\omega_{k}+P_{k}^{*}-\tilde{P}_{k}\\\
\tilde{P}_{k}&\approx P_{k}\\\
\tau_{f}\frac{d}{dt}\tilde{Q}_{k}&=Q_{k}-\tilde{Q}_{k}\\\
V_{m,k}&=V_{m,k}^{*}+R_{q}\big{(}Q_{k}^{*}-\tilde{Q}_{k}\big{)}\end{split}$
(1)

Here the dynamics are given as they are implemented in numerical integration.
In particular, five variables ($\omega_{k},\theta_{k},V_{m,k},Q_{k},P_{k}$)
are manipulated and three ultimate state variables
$\omega_{k},\theta_{k},V_{m,k}$ are passed further to the inner
voltage/current loops for subsequent physical implementation on the network
node. Power variables $Q_{k}$ and $P_{k}$ signify physical reactive/active
powers entering/leaving the node, while their counterparts $\tilde{Q}_{k}$ and
$\tilde{P}_{k}$ represent their averaged (low-pass filtered) measurements.
$\tau_{f}$ parameter in $Q_{k}$ relation is the low-pass time-constant. Note
that active power filter dynamics are neglected because the subsequent virtual
inertia dynamics usually are much longer in time ($M_{k}/D_{k}>>\tau_{f}$) so
that by the time power terms have any effect on $\omega_{k}$ dynamics, the
active power is already stabilized.

By eliminating the measurement dynamics one can reach to the system
representation of solely state variable dynamics (2):

$\begin{split}\frac{d}{dt}\theta_{k}&=\omega_{k}\\\
\tau_{f}\frac{d}{dt}V_{m,k}&=R_{q}\big{(}Q_{k}^{*}-Q_{k}\big{)}+\big{(}V_{m,k}^{*}-V_{m,k}\big{)}\\\
M_{k}\frac{d}{dt}\omega_{k}&=-D_{k}\omega_{k}+P_{k}^{*}-P_{k}\end{split}$ (2)

State dynamics of the form (2) can be referred to as reference nodal dynamics,
through which the regulator equivalence is established.

###  II-B Droop Control

Virtual nodal dynamics of the droop controller (as implemented numerically)
are given in (3):

$\begin{split}\frac{d}{dt}\theta_{k}&=\omega_{k}\\\
\omega_{k}&=R_{p}\big{(}P_{k}^{*}-\tilde{P}_{k}\big{)}\\\
\tau_{f}\frac{d}{dt}\tilde{P}_{k}&=P_{k}-\tilde{P}_{k}\\\
\tau_{f}\frac{d}{dt}\tilde{Q}_{k}&=Q_{k}-\tilde{Q}_{k}\\\
V_{m,k}&=V_{m,k}^{*}+R_{q}\big{(}Q_{k}^{*}-\tilde{Q}_{k}\big{)}\end{split}$
(3)

Here, instead of mimicking the synchronous machine relations, dq-frame angular
speed $\omega_{k}$ is made directly proportional to the estimated power
mismatch between the reference power and the nodal measurement. Therefore,
low-pass filter response on active power measurement can not be neglected.
From (3), one arrives at the state variable dynamics by eliminating
measurement variables:

$\begin{split}\frac{d}{dt}\theta_{k}&=\omega_{k}\\\
\tau_{f}\frac{d}{dt}V_{m,k}&=R_{q}\big{(}Q_{k}^{*}-Q_{k}\big{)}+\big{(}V_{m,k}^{*}-V_{m,k}\big{)}\\\
\underbrace{\frac{\tau_{f}}{R_{p}}}_{M_{dr,k}}\frac{d}{dt}\omega_{k}&=-\underbrace{\frac{1}{R_{p}}}_{D_{dr,k}}\omega_{k}+P_{k}^{*}-P_{k}\end{split}$
(4)

From (4) it can be seen that droop state dynamics have the same form as in VSM
case, thus proving equivalence. The corresponding droop inertia constant
$M_{dr,k}=\tau_{f}/R_{p}$ is the ratio of the filter time-constant and the
active power droop gain. The corresponding droop damping constant
$D_{dr,k}=1/R_{p}$ is the inverse of the active power droop gain. Similar
results were obtained in [6].

With the filter time-constant assumption $\tau_{f}<<M_{k}/D_{k}$ one can
conclude that for droop control the nodal angle/frequency dynamics are much
faster in converging to a steady-state compared to the VSM. On the other hand,
one would expect $M_{dr,k}<M_{k}$ so that the initial rate of change of
frequency (ROCOF) is greater in case of droop control, resulting in greater
transient frequency overshoots.

###  II-C Matching Control

Matching control is a relatively new alternative to the established
counterparts presented above. Some work on matching was presented in [7], [8].
The virtual/physical nodal dynamics of matching control are presented in (5):

$\begin{split}\frac{d}{dt}\theta_{k}&=\omega_{k}\\\
V_{DC,k}&=\tilde{V}_{DC,k}\\\
\tau_{f}\frac{d}{dt}\tilde{Q}_{k}&=Q_{k}-\tilde{Q}_{k}\\\
\omega_{k}&=K_{\theta}\big{(}\tilde{V}_{DC,k}-V_{DC,k}^{*}\big{)}\\\
V_{m,k}&=V_{m,k}^{*}+R_{q}\big{(}Q_{k}^{*}-\tilde{Q}_{k}\big{)}\\\
C_{DC,k}\frac{d}{dt}V_{DC,k}&=i_{DC,k}-\frac{1}{V_{DC,k}}P_{k}\end{split}$ (5)

Here, unlike techniques above, the frequency is directly related to the DC/AC
converter DC voltage mismatch rather than to the active power imbalance in
some form. Another differentiating feature is the fact that DC voltage
dynamics (last equation in (5)) is physical rather than virtual, meaning that
the dynamics describe physical capacitor charing/discharging rather than its
virtual counterpart. Alternatively, one could claim that the DC voltage
integration forward in time occurs ”physically” rather than numerically. As
before, the remaining dynamics are still virtual, implying the need for DC
voltage measurement. In (5), DC voltage measurement was approximated without
low-pas filter as in case of active power in VSM dynamics. In this case,
however, it is not due to subsequent dynamic time-scale separation, but due to
the fact that physical DC voltage dynamics already model averaged voltage,
neglecting the high-frequency switching ripple. AC voltage magnitude dynamics
remain unchanged.

Presented DC-side voltage dynamics have a DC current injection $i_{DC,k}$
coming from the preceding converter/power source and serving as a control
handle for DC voltage dynamics. Using this handle to establish a DC voltage
proportional control results in equivalent state nodal dynamics as in VSM and
droop as it will be shown in the following. DC voltage controller equations
are as follows (6):

$\begin{split}i_{DC,k}&=\tilde{i}_{DC,k}\\\
\tilde{i}_{DC,k}&=i_{DC,k}^{*}+K_{DC}\big{(}V_{DC,k}^{*}-\tilde{V}_{DC,k}\big{)}\end{split}$
(6)

Here $\tilde{i}_{DC,k}$ corresponds to the regulator output, while $i_{DC,k}$
is its physical manifestation. Solving together (5)-(6) for state dynamics as
before, one arrives at (7):

$\begin{split}\frac{d}{dt}\theta_{k}&=\omega_{k}\\\
\tau_{f}\frac{d}{dt}V_{m,k}&=R_{q}\big{(}Q_{k}^{*}-Q_{k}\big{)}+\big{(}V_{m,k}^{*}-V_{m,k}\big{)}\\\
\underbrace{\frac{C_{DC,k}V_{DC,k}^{*}}{K_{\theta}}}_{M_{DC,k}}\frac{d}{dt}\omega_{k}&=-\underbrace{\frac{K_{DC}V_{DC,k}^{*}}{K_{\theta}}}_{D_{DC,k}}\omega_{k}+...\\\
&...+\underbrace{V_{DC,k}^{*}i_{DC,k}^{*}}_{P_{k}^{*}}-\underbrace{\frac{V_{DC,k}^{*}}{V_{DC,k}}}_{\approx
1}P_{k}\end{split}$ (7)

Again, from (7) one can observe equivalence to the VSM/droop state dynamics
with $M_{DC,k}$ and $D_{DC,k}$ depending on DC-link physical parameters
$C_{DC,k},V_{DC,k}^{*}$ as well as control gains $K_{DC},K_{\theta}$.

Using the control gains one could shape the equivalent inertia and damping
parameters to achieve the desired nodal transient performance. For instance,
to limit the ROCOF one could reduce the $K_{\theta}$ gain, thus increasing the
equivalent inertia $M_{DC,k}$. To reduce the steady state convergence time one
could increase $K_{DC}$ gain. It would reduce the transient time-constant
$\tau_{DC,k}=M_{DC,k}/D_{DC,k}=C_{DC,k}/K_{DC}$ and lead to faster convergence
of nodal states.

Next section provides more insights on desired nodal parameters and gains. As
of now, the equivalent inertias and damping parameters of three grid-forming
controllers are summarized in Table I.

TABLE I: Equivalent Inertia and Damping Parameters Parameter | Control Technique
---|---
$\big{(}M_{k},D_{k}\big{)}$ | VSM | Droop | Matching
$M_{k}$ | $M_{k}$ | $\frac{\tau_{f}}{R_{p}}$ | $\frac{C_{DC,k}V_{DC,k}^{*}}{K_{\theta}}$
$D_{k}$ | $D_{k}$ | $\frac{1}{R_{p}}$ | $\frac{K_{DC}V_{DC,k}^{*}}{K_{\theta}}$

##  III Network Dynamics of Grid-Forming Controls

Having established the equivalence among the nodal dynamics of the grid-
forming controllers under consideration, this section addresses the network
dynamics of the graph with power converters on its nodes and transmission
lines across its edges. Transmission line AC power flows are summarized in the
following equation (8):

$\begin{split}P_{k}&=\sum_{l=1}^{n}||V_{k}||\cdot||V_{l}||\cdot\big{[}G_{kl}cos(\theta_{k}-\theta_{l})+B_{kl}sin(\theta_{k}-\theta_{l})\big{]}\\\
Q_{k}&=\sum_{l=1}^{n}||V_{k}||\cdot||V_{l}||\cdot\big{[}G_{kl}sin(\theta_{k}-\theta_{l})+B_{kl}cos(\theta_{k}-\theta_{l})\big{]}\end{split}$
(8)

Dealing with these highly nonlinear functions of nodal angle differences poses
a significant challenge. To simplify the theoretical analysis, AC power flow
is transformed to so called DC power flow by linearizing the AC equations
around the nominal trivial operating point
$P_{k}^{*}=P_{k}=Q_{k}^{*}=Q_{k}=0,\forall k$. Resulting linearized power
flows are as follows:

$\begin{split}P_{k}&=\sum_{l=1}^{n}B_{kl}\big{(}\theta_{k}-\theta_{l}\big{)}\\\
Q_{k}&=\sum_{l=1}^{n}B_{kl}\big{(}V_{m,k}-V_{m,l}\big{)}\end{split}$ (9)

Further simplifying assumption is that all converters are tuned identically,
meaning the same $M_{k}=m,D_{k}=d,V_{m,k}^{*}=V_{m}^{*}$ for $\forall k$.
Through linearization and vectorization of nodal equation (2) for all nodes
$k\in[1,..,n]$ one arrives at linearized VSM network dynamics in dq-frame
(10):

$\begin{split}\frac{d}{dt}\theta_{\delta}&=\omega_{\delta}\in\mathbb{R}^{n\times
1}\\\
\frac{d}{dt}\omega_{\delta}&=-\frac{1}{m}L_{B}\theta_{\delta}-\frac{d}{m}\omega_{\delta}\in\mathbb{R}^{n\times
1}\\\
\frac{d}{dt}V_{m\delta}&=-\frac{1}{\tau_{f}}\big{[}L_{B}R_{q}+I_{n}\big{]}V_{m\delta}\in\mathbb{R}^{n\times
1}\end{split}$ (10)

Representing the first two dynamics as a system (11):

$\frac{d}{dt}\begin{bmatrix}\theta_{\delta}\\\
\omega_{\delta}\end{bmatrix}=\underbrace{\begin{bmatrix}0_{n\times n}&I_{n}\\\
-\frac{L_{B}}{m}&-\frac{d}{m}I_{n}\end{bmatrix}}_{\text{$\mathcal{L}\in\mathbb{R}^{2n\times
2n}$}}\begin{bmatrix}\theta_{\delta}\\\
\omega_{\delta}\end{bmatrix}\in\mathbb{R}^{2n\times 1}$ (11)

By $L_{B}\in\mathbb{R}^{n\times 1}$ one denotes the Laplacian matrix of
susceptances $B$ of the network (graph). By $m$ and $d$ one denotes inertia
and damping constants of any single power converter an any node. Using the
equivalence from the last section, subsequent discussion is valid for all
three control techniques under consideration depending on the values of $m$
and $d$ selected.

Analyzing the system (11) one can claim the following:

  * •

Laplacian matrix $L_{B}$ has $n$ eigenvalues $\lambda(L)$. One of them is
strictly at 0, others are real and have strictly positive real parts. The more
the connectivity of the corresponding graph of $L_{B}$ the greater the non-
zero eigenvalues are.

  * •

Eigenvalues of the so called ”larger” Laplacian matrix $\mathcal{L}$ depend on
eigenvalues $\lambda(L)$ as follows (12):

$\eta=-\frac{d}{2m}\pm\frac{1}{2}\sqrt{\frac{d^{2}}{m^{2}}-\frac{4\lambda}{m}}$
(12)

So there is one eigenvalue at $\eta=0$, one at $\eta=-d/m$ and $2n-2$
potentially complex eigenvalues with strictly negative real parts.

  * •

Using appropriate matrix similarity transformation $T\in\mathbb{R}^{2n\times
2n}$, the dynamics in (11) can be expressed in the form of average angle/speed
and respective difference vectors
$[\theta_{AVG};\omega_{AVG};\Delta\theta_{\delta};\Delta\omega_{\delta}]$.
Next, it can be shown that $\eta=0$ eigenvalue corresponds to the
$d\theta_{AVG}/dt=\omega_{AVG}$ dynamics. The remaining state dynamics do not
depend on $\theta_{AVG},\omega_{AVG}$, so the corresponding
$L^{\prime}\in\mathbb{R}^{2n-2\times 2n-2}$ is Hurwitz, implying asymptotic
stability of the states $[\Delta\theta_{\delta};\Delta\omega_{\delta}]$ around
$[0;0]$

  * •

It may be shown that by adding the power disturbance of the form
$P_{d}=[P_{d,AVG};\Delta P_{d}]$ angle difference vector
$\Delta\theta_{\delta}$ and average speed $\omega_{AVG}$ converge to some non-
zero disturbance-dependent constants. Average angle $\theta_{AVG}$ continues
integrating as a result of non-zero $\omega_{AVG}$ and corresponding
eigenvalue $\eta=0$.

  * •

Second largest eigenvalue of $\mathcal{L}$, $\eta_{2}$ determines the rate of
convergence of the states to steady state. From (12) the smaller (more
negative) the $\eta_{2}$ the more connected the graph, the faster the
convergence and potentially the more oscillatory the transient. Indeed,
smaller $\eta_{2}$ implies greater $\lambda$’s - connected graph. Also smaller
$\eta_{2}$ means transient exponential decaying at faster rate. Finally,
greater $\lambda$’s can lead to imaginary components in $\eta$’s (negative
values under square root), resulting in oscillatory (underdamped) transient
decay of the states.

From these it can be deduced that in order to have the desired network-level
dynamics (at least within the idealized setting at hand), one needs to have a
large $d/m$ ratio, such that the eigenvalues $\eta$’s are far into the left
half-plane and strictly real even for highly connected graph (high
$\lambda$’s). This, in turn, means high damping constant $d$ and low inertia
constant $m$, resembling droop parameters.

On the other hand, the ROCOF has to be kept within specified standards to
comply with practical limits of protection devices. Thus inertia parameters
$m$ cannot be too low. Overall, the desired parameter combination is low
enough $m$ (to limit ROCOF) and a high $d$ (to stabilize and accelerate the
transient even for highly connected networks). Table I can be used for
translating these requirements to the gains of respective controllers.

Last but not least, despite the fact that equivalence was established among
the three grid-forming controllers at hand, there are practical advantages
associated with matching control that make it an attractive option for a
further research. First, matching control does not assume a constant stiff DC
source behind the inverter. Second, instead of controlling active power
directly, matching control regulates the DC voltage and attaches $\omega$ to
it. As it was observed, the net result is the same nodal grid-forming
($\theta,\omega,V$-forming) operation of a power converter. However, matching
does a step in the direction of right causality, assuming that some power
source behind the DC-link is the main generator of energy, not the inverter
itself as in case of VSM or droop control.

##  IV Conclusion and Outlook

To sum up, this work addressed the equivalence of three grid-forming control
techniques (VSM, droop and matching) using theoretical analysis of a
simplified nodal converter and network system. For nodal analysis the AC
voltage/current controllers are assumed to be ideal. Also in case of VSM/droop
the DC voltage was assumed to be stiff. For network analysis, all node
converters are assumed to be identical and the AC power flow was linearized
around trivial operating point. These simplifications allow the theoretical
study of the grid-forming controllers at hand. Through this analysis the
equivalence of the techniques was shown both in nodal and network-level
dynamics. Droop and matching counterparts of the VSM dynamic parameters were
presented. Network-level analysis showed how the VSM-equivalent parameters
should be selected to achieve network-level objectives of a stable, fast and
non-oscillatory transients while respecting the nodal constraints, such as
rate of change of frequency (ROCOF). Overall, with established equivalence of
the techniques, one would suggest to use matching control as it incorporates
the DC voltage dynamics and makes a step towards a right power flow causality.

As an outlook, one foresees the work in the direction of matching control.
Establishing fully correct causality between grid and RES, current limiting
during transients, network planning issues - to name just a few challenges yet
to be addressed.

## References

  * [1] B. Kroposki et al., ”Achieving a 100% Renewable Grid: Operating Electric Power Systems with Extremely High Levels of Variable Renewable Energy,” in IEEE Power and Energy Magazine, vol. 15, no. 2, pp. 61-73, March-April 2017, doi: 10.1109/MPE.2016.2637122.
  * [2] Q. Peng, Q. Jiang, Y. Yang, T. Liu, H. Wang and F. Blaabjerg, ”On the Stability of Power Electronics-Dominated Systems: Challenges and Potential Solutions,” in IEEE Transactions on Industry Applications, vol. 55, no. 6, pp. 7657-7670, Nov.-Dec. 2019, doi: 10.1109/TIA.2019.2936788.
  * [3] J. Rocabert, A. Luna, F. Blaabjerg and P. Rodríguez, ”Control of Power Converters in AC Microgrids,” in IEEE Transactions on Power Electronics, vol. 27, no. 11, pp. 4734-4749, Nov. 2012, doi: 10.1109/TPEL.2012.2199334.
  * [4] B. K. Poolla, D. Groß and F. Dörfler, ”Placement and Implementation of Grid-Forming and Grid-Following Virtual Inertia and Fast Frequency Response,” in IEEE Transactions on Power Systems, vol. 34, no. 4, pp. 3035-3046, July 2019, doi: 10.1109/TPWRS.2019.2892290.
  * [5] A. Tayyebi, D. Groß, A. Anta, F. Kupzog and F. Dörfler, ”Frequency Stability of Synchronous Machines and Grid-Forming Power Converters,” in IEEE Journal of Emerging and Selected Topics in Power Electronics, vol. 8, no. 2, pp. 1004-1018, June 2020, doi: 10.1109/JESTPE.2020.2966524.
  * [6] S. D’Arco and J. A. Suul, ”Equivalence of Virtual Synchronous Machines and Frequency-Droops for Converter-Based MicroGrids,” in IEEE Transactions on Smart Grid, vol. 5, no. 1, pp. 394-395, Jan. 2014, doi: 10.1109/TSG.2013.2288000.
  * [7] S. Curi, D. Groß and F. Dörfler, ”Control of low-inertia power grids: A model reduction approach,” 2017 IEEE 56th Annual Conference on Decision and Control (CDC), Melbourne, VIC, 2017, pp. 5708-5713, doi: 10.1109/CDC.2017.8264521.
  * [8] C. Arghir and F. Dörfler, ”The Electronic Realization of Synchronous Machines: Model Matching, Angle Tracking, and Energy Shaping Techniques,” in IEEE Transactions on Power Electronics, vol. 35, no. 4, pp. 4398-4410, April 2020, doi: 10.1109/TPEL.2019.2939710.