11institutetext:  Leiden Observatory, Leiden University, P.O. Box 9513, 2300
RA Leiden, The Netherlands
11email<EMAIL_ADDRESS>22institutetext: ASTRON Netherlands
Institute for Radio Astronomy, Oude Hoogeveensedijk 4, 7991 PD Dwingeloo, The
Netherlands 33institutetext: School of Physical Sciences and Centre for
Astrophysics & Relativity, Dublin City University, Glasnevin, D09 W6Y4,
Ireland 44institutetext: Centre for Astrophysics Research, Department of
Physics, Astronomy and Mathematics, University of Hertfordshire, College Lane,
Hatfield AL10 9AB, UK 55institutetext: GEPI & USN, Observatoire de Paris,
Université PSL, CNRS, 5 Place Jules Janssen, 92190 Meudon, France
66institutetext: Department of Physics & Electronics, Rhodes University, PO
Box 94, Grahamstown, 6140, South Africa 77institutetext: Anton Pannekoek
Institute for Astronomy, University of Amsterdam, Science Park 904, 1098 XH
Amsterdam, The Netherlands

# Possible discovery of Calvera’s supernova remnant

M. Arias 11 A. Botteon 11 C. G. Bassa 22 S. van der Jagt 11 R. J. van Weeren
11 S. P. O’Sullivan 33 Q. Bosschaart 11 R. S. Dullaart 11 M. J. Hardcastle 44
J. W. T. Hessels 2255 T. Shimwell 22 M. M. Slob 11 J. A. Sturm  C. Tasse 5566
N. C. M. A. Theijssen 11 J. Vink 77

(Received ; Accepted )

We report the discovery of a ring of low surface brightness radio emission
around the Calvera pulsar, a high Galactic latitude, isolated neutron star, in
the LOFAR Two-metre Sky Survey (LoTSS). It is centered at
$\alpha=14\mathrm{h}11\mathrm{m}12\aas@@fstack{s}6$, $\delta=+79\degr
23\arcmin 15\arcsec$, has inner and outer radii of $14\aas@@fstack{\prime}2$
and $28\aas@@fstack{\prime}4$, and an integrated flux density at 144 MHz of
$1.08\pm 0.15$ Jy. The ring center is offset by $4\aas@@fstack{\prime}9$ from
the location of the Calvera pulsar. H$\alpha$ observations with the Isaac
Newton Telescope show no coincident optical emission, but do show a small
($\sim 20$″) optical structure internal to the ring. We consider three
possible interpretations for the ring: that it is an H II region, a supernova
remnant (SNR), or an Odd Radio Circle (ORC). The positional coincidence of the
ring, the pulsar, and an X-ray-emitting non-equilibrium ionisation plasma
previously detected, lead us to prefer the SNR interpretation. If the source
is indeed a SNR and its association with the Calvera pulsar is confirmed, then
Calvera’s SNR, or G118.4+37.0, will be one of few SNRs in the Galactic halo.

###### Key Words.:

Radio continuum: general – ISM: supernova remnants – ISM: HII regions – Stars:
neutron – pulsars: individual (1RXS J141256.0+792204)

##  1 Introduction

The LOFAR Two-metre Sky Survey (LoTSS, Shimwell et al., 2017) is an ongoing
survey of the Northern sky at $120-168$ MHz, with $6\arcsec$ resolution and
$\lesssim 100$ $\upmu$Jy beam-1 sensitivity111The LoTSS data are taken with
the full configuration of the LOFAR array, including international stations,
but only the Dutch configuration of the array is processed. The data are
available for reprocessing to $0.3\arcsec$ resolution.. There have been two
LoTSS data releases (Shimwell et al., 2019, 2022) that have enabled a wide
variety of scientific results.

While processing proprietary data from LoTSS, we have discovered a ring of
diffuse radio emission encompassing the celestial position of the isolated
neutron star Calvera (1RXS J141256.0+792204). It is a degree-scale circular
structure at high Galactic latitude ($l=37\degr$), that may or may not be
associated with the Calvera neutron star.

1RXS J141256.0+792204 is an isolated neutron that so far defies classification
into the established types of neutron stars. Initially identified as a
candidate thermally emitting isolated neutron star by Rutledge et al. (2008),
selected on its soft X-ray spectrum and high X-ray to optical flux ratio from
the ROSAT All-Sky Survey (Voges et al., 1999), it was nicknamed “Calvera” due
to its similarities with the 7 nearby X-ray dim isolated neutron stars
(XDINS), colloquially named the “Magnificent Seven” (Haberl, 2007). However,
its high high black body temperature set it apart from the XDINS, a
distinction further confirmed by the discovery of X-ray pulsations at a 59 ms
spin period (Zane et al., 2011), much shorter than the several second spin
periods of the “Magnificent Seven”, and by a distance estimate of 3.3 kpc
(Mereghetti et al., 2021), much larger than the few hundred parsec distances
of the other XDINS. X-ray timing yields dipolar spindown properties
($\dot{E}=6.1\times 10^{35}$ erg s-1, $\tau_{\mathrm{c}}=2.9\times 10^{5}$ yr,
$B=4.4\times 10^{11}$ G) that are more in line with those of young radio
pulsars (Halpern et al., 2013; Bogdanov et al., 2019; Mereghetti et al.,
2021), disfavouring the classification as a nearby millisecond pulsar or a
central compact object (CCO). On the other hand, the absence of radio
pulsations (Hessels et al., 2007; Zane et al., 2011) or gamma-ray pulsations
(Halpern, 2011; Halpern et al., 2013; Mereghetti et al., 2021) are somewhat
inconsistent with those of young radio pulsars.

Here, we report the detection of a diffuse radio ring around the Calvera
pulsar, and discuss three possible scenarios: that it is an H II region, that
it is a supernova remnant (SNR), possibly associated with Calvera, or that it
is an Odd Radio Circle (ORC).

Figure 1: a) Radio emission from the diffuse ring at 144 MHz in the LOFAR HBA
LoTSS imaging mosaic ($1\aas@@fstack{\circ}25\times 1\aas@@fstack{\circ}25$
field-of-view, $20\arcsec$ resolution). The location of Calvera is indicated
with a circle, with the arrow representing the Halpern & Gotthelf (2015)
proper motion from 10 kyr in the past to 10 kyr into the future. The dashed
lines indicate the propagated $1\sigma$ positional uncertainty due to proper
motion. The cross sign indicates the estimated center of the diffuse ring,
with the dotted circles representing the inner and outer radii. The dashed
circle ($4\arcmin$ radius) denotes the location of the diffuse X-ray emission
seen by Zane et al. (2011) from 1RXS J140818.1+792113, with the square
indicating the location of H$\alpha$ emission. b) The diffuse ring in the
Westerbork Synthesis Radio Telescope WENSS survey imaging at 325 MHz and c)
the Very Large Array NVSS survey at 1.4 GHz. Both panels have a
$1\aas@@fstack{\circ}25\times 1\aas@@fstack{\circ}25$ field-of-view and
identical pointing and orientation to that of panel a. d) The 0.1-2.4 keV
X-ray image from the ROSAT All-Sky Survey (Voges et al., 1999), with the ring
center and inner outer radii of the radio ring denoted by the cross and dotted
circles. The dashed circle encircles the emission from 1RXS J140818.1+792113.
e) A $2\arcmin\times 2\arcmin$ section of the $6\arcsec$ resolution LoTSS
imaging mosaic at the location of Calvera, indicated with the dotted circle
($10\arcsec$ radius).

##  2 Observations and analysis

###  2.1 Radio continuum imaging

Visual inspection of radio continuum images obtained as part of the ongoing
LOFAR Two-metre Sky Survey (LoTSS; Shimwell et al. 2017) led to the discovery
of a $\sim 1\degr$ diameter diffuse ring-like structure at
$l=118\aas@@fstack{\circ}4$, $b=+37\aas@@fstack{\circ}0$ in the 8 hr LoTSS
pointing P207+77 (observed on 2019 September 19). The presence of the ring was
confirmed in two nearby LoTSS pointings, P207+80 and P222+80, observed on 2022
January 29 and 2022 February 6, respectively, both with 8 hr exposures. The
LoTSS pointings were offset from the center of the ring by
$2\aas@@fstack{\circ}20$ (P207+77), $1\aas@@fstack{\circ}04$ (P207+80) and
$1\aas@@fstack{\circ}93$ (P222+80). These pointings were processed by the
LOFAR Surveys Key Science Project team, employing the ddf-
pipeline222https://github.com/mhardcastle/ddf-pipeline developed to reduce
LoTSS observations (Shimwell et al., 2017, 2019, 2022), which delivers
direction-dependent corrected mosaics at $6\arcsec$ and $20\arcsec$ resolution
(see Tasse et al., 2021, for more details). The imaging is done using the
multiscale cleaning algorithm described in Offringa & Smirnov (2017).

Since the edges of the ring are not sharp, we estimate its center, radius and
width by modelling it as a circularly symmetric radially offset Gaussian.
Starting with a $1\aas@@fstack{\circ}25\times 1\aas@@fstack{\circ}25$ section
centered on the ring from the $20\arcsec$ resolution LoTSS image, all pixels
above the 98% percentile (0.7 mJy beam-1) were masked. The masked shapes were
expanded by one pixel through a single binary dilation operation to mask the
fainter wings of point sources. In the resulting image, 3.1% of the pixels
were masked. The unmasked pixels were used in a least-squares minimisation of
the circularly symmetric Gaussian, fitting for the Gaussian width, the radial
offset around the center, as well as the Gaussian amplitude above some
constant background.

To measure the integrated flux density of the ring, we further processed the
LoTSS data using the “extraction and selfcal” method described in van Weeren
et al. (2021). This methods allows one to correct for residual calibration
artifacts while permitting fast and flexible reimaging of a smaller portion of
the dataset. We focused on a squared region of 1.1 deg2 centered on the ring
and obtained a new image with WSClean (Offringa et al., 2014) at 30″
resolution. In this image, obtained using a minimum uv-cut of $60\lambda$
(corresponding to an angular size of $\sim 57\arcmin$), discrete sources were
subtracted from the visibilities using the clean components created from a
previous imaging run at higher resolution (created with an minimum uv-cut of
$1000\lambda$ to filter out the ring).

###  2.2 Radio Polarisation

The LoTSS data products include linear polarisation Stokes Q and Stokes U
cubes with a spatial resolution of 20″. To search for linearly polarized
emission, we used the rotation measure synthesis technique (Brentjens & de
Bruyn, 2005) on the 3′ resolution Q and U channel images with
pyrmsynth_lite333https://github.com/sabourke/pyrmsynth_lite. There were 480 Q
and U images from 120 MHz to 168 MHz. We searched the band-averaged polarized
intensity images across a Faraday depth range of $\pm 100$ rad m-2 with a
sampling of 0.5 rad m-2 and found no evidence for polarized emission
associated with the ring, at a level of 300 $\mu$Jy beam-1. This is a higher
noise than typical for LoTSS images, but is due to the fact that the ring is
located at the edge of the observing field. More details on the rotation
measure synthesis method developed for LoTSS data products can be found in
O’Sullivan et al. (in prep).

###  2.3 Beamformed radio observations

A high-time resolution Stokes I observation (ObsID: L257877) of Calvera was
obtained on 2015 January 16 using the inner 22 LOFAR HBA core stations in
beamforming mode, providing a $3\aas@@fstack{\prime}5$ tied-array beam pointed
at the position of Calvera. The observation was 3 hr in length with 78.125 MHz
of bandwidth, centred at 149 MHz, and with 163.84 $\upmu$s time resolution and
12.2 kHz frequency resolution. For this observational setup, smearing due to
dispersion reaches 1 ms at a dispersion measure (DM) of about 30 pc cm-3 and 3
ms at $\mathrm{DM}=80$ pc cm-3. Models for the Galactic electron distribution
predict maximum dispersion measures (DMs) of 42.2 pc cm-3 at 10.1 kpc (NE2001;
Cordes & Lazio 2002) and 35.5 pc cm-3 at 15.2 kpc (YMW16, Yao et al. 2017)

The observation was folded with the Halpern & Gotthelf (2015) X-ray timing
ephemeris using dspsr (van Straten & Bailes, 2011) to 10 s sub-integrations.
Though the observation epoch falls outside of the timing span of both the
Halpern & Gotthelf (2015) and Mereghetti et al. (2021) timing ephemerides,
extrapolating these to the observation epoch yields consistent folding periods
($\Delta P=-17.6\pm 6.9$ ns). Using psrchive (Hotan et al., 2004) tools, radio
frequency interference (RFI) was masked and the folded observation was
subsequently dedispersed to 401 trial dispersion measures (DMs) between 0 and
80 pc cm-3 at 0.2 pc cm-3 steps. Each DM trial was averaged to 128 pulse phase
bins, 0.78 MHz frequency channels and 100 s sub-integrations and searched for
periodicities within 280 ns around the folding period and DM offsets of
$-0.23$ to $0.23$ pc cm-3 with pdmp. No periodicities with a significance over
$7\sigma$ were found. The flux calibration method by Kondratiev et al. (2016)
estimates sensitivity limits of 0.25 mJy for an assumed fractional pulse width
of 10% (0.8 mJy for 50%) at $7\sigma$.

The observation was also searched for emission from single pulses using tools
from the presto software suite (Ransom, 2001). After RFI masking, dedispersed
timeseries were created for DMs up to 80 pc cm-3, with steps of 0.002 pc cm-2
up to 10 pc cm-3 at the native time resolution, 0.005 pc cm-3 from 10 to 22.5
pc cm-3 for a downsampled time resolution of 0.655 ms, and 0.02 pc cm-3 at
1.31 ms. Single pulse candidates with DMs below 0.2 pc cm-3 were discarded due
to pollution by zero DM RFI, while all other events above a significance of
$7\sigma$ were visually inspected. Unfortunately, no astrophysical pulses were
detected as all remaining events could be attributed to time ranges affected
by RFI. We estimate $7\sigma$ fluence limits of 2.6 Jy ms for 1 ms wide
signals.

Figure 2: The detection of extended H$\alpha$ emission in the INT/WFC
observations in the narrowband H$\alpha$ filter (panel a), the broadband
$r$-band filter (panel b) and the continuum subtracted H$\alpha-r$-band (panel
c). Each panel shows a $4\arcmin\times 4\arcmin$ field-of-view, with the
dashed circle representing extent of the diffuse X-ray emission from 1RXS
J140818.1+792113 (see Fig. 1). Radio contours from the $6\arcsec$ LOFAR mosaic
are overlaid in panel a.

###  2.4 Optical Observations

We also observed the region around Calvera and a section of the radio ring
with the the Wide Field Camera (WFC) on the 2.5m Isaac Newton Telescope (INT)
in La Palma. Our initial aims were to search for optical forbidden line
emission from the radio shell that would indicate radiative SNR shocks, and
the presence of a bow shock near Calvera that could help us constrain its
direction of motion. The region was observed for three nights (2022 March 19,
20, and 23) in the H$\alpha$ and $r$-band filters.

In order to reduce the optical observing data we used the automated
astronomical data reduction pipeline from THELI (Erben et al., 2005; Schirmer,
2013). THELI provided the astrometry and relative photometry to combine the
observations of a single observation night. For every single observation THELI
subtracted the sky with a model. To model the sky THELI uses SExtractor
(Bertin & Arnouts, 1996) to remove all sources, dynamically fills the gaps in
the field and convolves the resulting field with a Gaussian kernel of 256
pixels. Finally we used THELI to coadd the observations of a single
observation night and perform a flux calibration with known field sources from
Pan-STARRS (Chambers et al., 2016a). For the $r$-band and H$\alpha$ we also
combined the observations for different nights. The $r$-band and H$\alpha$
images were photometrically calibrated against the PanSTARSS $r$-band
magnitudes (Chambers et al., 2016b) using THELI.

No emission associated with either the radio ring or a bow shock was detected
in H$\alpha$ (174 minutes of exposure time). We estimate a $5\sigma$ limit on
H$\alpha$ surface density flux of $1.1\times 10^{-16}$ erg cm-2 s-1 arcsec-2.

We did detect, though, a smudge of H$\alpha$ emission centered at
$\alpha_{\mathrm{J2000}}=14^{\mathrm{h}}08^{\mathrm{m}}28\aas@@fstack{s}3$,
$\delta_{\mathrm{J2000}}=+79\degr 22\arcmin 54\aas@@fstack{\prime\prime}0$.
The extended emission is faint, with an H$\alpha$ surface brightness of 24.0
mag arcsec-2, somewhat V-shaped, and at its longest it is $23\arcsec$ across.
Figure 2 shows the H$\alpha$ and $r$-band images, as well as the continuum
subtracted H$\alpha$ image for a region surrounding the smudge. There is no
clear radio emission associated with the H$\alpha$ smudge, although there is a
point source that almost overlaps with the smudge, and so it is hard to tell
if faint extended emission is also present. We examined the four WISE bands
for emission in this region and found no clear IR counterpart. The smudge is
unlikely to be of extra-galactic origin for several reasons: its angular
scale, the fact that it does not overlap in full with a galaxy, and the fact
that it falls within the redshift range of the filter ($z<0.007$, for an
H$\alpha$ filter bandpass of 95
Å444https://www.ing.iac.es/astronomy/instruments/wfc/).

##  3 Results

###  3.1 Detection of a radio ring

We report the detection of radio emission from a degree-scale diffuse ring in
LoTSS radio continuum imaging at 144 MHz. Figure 1a shows the diffuse ring in
the LoTSS radio continuum mosaic. The ring is located at
$\alpha=14^{\mathrm{h}}11^{\mathrm{m}}12\aas@@fstack{s}6$, $\delta=+79\degr
23\arcmin 15\arcsec$, corresponding to a Galactic position of
$l=118\aas@@fstack{\circ}41$, $b=37\aas@@fstack{\circ}03$. The ring has a FWHM
of $8\aas@@fstack{\prime}5$ and we conservatively estimate an the inner
diameter of $28\aas@@fstack{\prime}4$ and an outer diameter of
$56\aas@@fstack{\prime}8$ (see Fig. 1). These correspond to a physical inner
and outer diameter of 8.2 pc and 16.5 pc for a distance of 1 kpc. The ring
does not seem to have a sharp inner or outer edge, and its brightness is
relatively constant, both radially and in position angle, except for a
brighter, straight filament in the North-West. After subtracting the point
sources from the LoTSS mosaic, we measure an integrated flux density of
$1.08\pm 0.15$ Jy at 144 MHz (corresponding to the area between the dotted
concentric circles in Fig. 1). The uncertainty includes the statistical error,
and the systematic error on the flux scale (10% of the measured flux density
of the ring, following LoTSS Data Release 2, Shimwell et al., 2022).

Survey |  | Frequency
---
(MHz)
| Flux density
---
(Jy)
| Error in flux
---
density (Jy)
LoTSS | 144 | 1.08 | 0.15
WENSS | 325 | 0.63 | 0.21
NVSS | 1420 | 0.22 | 0.04
Table 1: Integrated flux density of the ring centered at
$\alpha=14\mathrm{h}11\mathrm{m}12\aas@@fstack{s}6$, $\delta=+79\degr
23\arcmin 15\arcsec$, with inner and outer radii of $14\aas@@fstack{\prime}2$
and $28\aas@@fstack{\prime}4$ at various frequencies.

Knowing the location and extent of the ring, we searched for its presence in
other radio continuum surveys. Radio emission of the diffuse ring is detected
at low significance in the Westerbork Northern Sky Survey (WENSS, Rengelink et
al., 1997) at 325 MHz (Fig. 1b), and at even lower significance in the NRAO
VLA Sky Survey (NVSS, Condon et al., 1998) at 1.4 GHz (Fig. 1c). The ring is
not detected in the NRAO VLA Low-frequency Sky Survey (VLSS; Cohen et al.
2007; Lane et al. 2014) at 74 MHz (0.7 Jy beam-1 sensitivity) and the TIFR
GMRT Sky Survey (TGSS; Intema et al. 2017) at 150 MHz (5 mJy beam-1
sensitivity). The integrated flux densities of the ring in the surveys where
it is detected are summarised in table 1, and correspond to a best-fit
spectral index value of $-0.71\pm 0.09$ (for power-law radio spectra with
observing frequency $\nu$ of the form $S_{\nu}\propto\nu^{\alpha}$).

The full area encompassed by the radio ring has only been observed in the
X-rays as part of the ROSAT All-Sky Survey (Voges et al., 1999)555It has
likely also been observed by eROSITA (Sunyaev et al., 2021), but the
observations are not public as of this writing., and no X-ray emission
coincident with the radio emission is detected (see Fig. 1d). Two ROSAT X-ray
sources from the bright source catalog (Voges et al., 1999) are coincident
with the extent of the ring, the Calvera X-ray pulsar (1RXS J141256.0+792204)
and 1RXS J140818.1+792113. Zane et al. (2011) find that the latter source is
extended and has a soft spectrum in XMM-Newton observations covering the inner
region of the ring, and that its X-ray spectrum best modelled as a non-
equilibrium ionization (NEI) plasma.

The Calvera X-ray pulsar is offset by $4\aas@@fstack{\prime}9$ from the center
of the radio ring, well within its inner radius. No other radio, X-ray or
gamma-ray pulsar from the ATNF Pulsar Catalogue (Manchester et al. 2005,
version 1.67) is located within the ring. Another 8 radio pulsars are within a
$10\degr$ radius, with the closest at an offset of $2\aas@@fstack{\circ}85$.

###  3.2 Non-detection of Calvera in LOFAR observations

No radio emission is detected in the LoTSS radio continuum mosaic at the
location of Calvera (Fig. 1e), either as a point source from pulsed or
unpulsed emission, or as extended emission from a pulsar wind nebula (PWN).
The rms noise measured at the location of Calvera in the 6″ resolution map is
75 $\mu$Jy beam-1; therefore, we estimate a $3\sigma$ flux density upper limit
of 225 $\mu$Jy for Calvera’s radio continuum emission at 144 MHz.

The beamformed LOFAR observations also fail to detect pulsed radio emission at
the known spin period of Calvera, or as single pulses down to flux density and
fluence limits of 0.25 mJy and 2.6 Jy ms at 149 MHz, respectively. These
limits add to previous pulsed flux density upper limits of 4 mJy at 385 MHz
and 0.3 mJy at 1380 MHz by Hessels et al. (2007) using the Westerbork
Synthesis Radio Telescope and 0.05 mJy at 1360 MHz with the Effelsberg
telescope (Zane et al., 2011). These observations confirm that Calvera is a
radio quiet pulsar.

###  3.3 Statistical modelling of the probability of chance alignment

Although chance alignments between pulsars and SNRs are frequent in the
Galactic plane (Gaensler & Johnston, 1995), at such high Galactic latitudes
both pulsars and remnants are rare. We performed some Monte Carlo simulations
so as to quantify the likelihood of a chance alignment.

In order to find the probability that a pulsar falls within a fixed area with
the coordinates and diameter of our radio ring (which we take to be one
degree), we simulate many iterations of pulsar populations and count the
number of times that a pulsar falls within the radio circle. Each sample
pulsar is generated with a set of Galactocentric coordinates $(r,\theta,h)$,
where $r$ is the Galactocentric radius, $\theta$ is the polar angle, and $h$
is the height above or below the Galactic plane. We sample $r$ according to
the distribution given by Yusifov & Küçük (2004) in Galactocentric radius, and
$\theta$ by drawing from a uniform distribution between 0 and $2\pi$. To
sample $h$, we tried using exponential distributions with scale heights of 100
pc, 330 pc and 350 pc (Faucher-Giguère & Kaspi, 2006; Lorimer et al., 2006;
Mdzinarishvili & Melikidze, 2004) and found that the simulated sample, once we
converted to heliocentric coordinates, was very compact in Galactic latitudes.
In runs of 10,000 Monte Carlo iterations, the rate at which a pulsar fell
within the area of the radio ring was consistently less than $5\times
10^{-4}$. This distribution produced on average 1 pulsar within a 10$\degr$
radius centered at Calvera, whereas we observe 9 in such a region, which
suggests that the scale height distributions in the above references do not
describe well this local region of the Galaxy.

Given that the simulated populations sampled from the exponential
distributions seemed to generate pulsars whose Galactic latitudes were too
low, we also simulated the distribution in $h$ empirically. We used the ATNF
Pulsar Catalogue v1.67 (Manchester et al., 2005), removing pulsars in globular
clusters or in the Magellanic Clouds, in order to histogram the distribution
in heights above the Galactic plane (using the dispersion measure distances).
Then, we sampled from the midpoints of the histogram bins, adding some noise
so as to not sample in discrete heights above or bellow the plane. Naturally,
since we sampled from the ATNF Pulsar Catalogue v1.67 empirical distribution,
we recover the observed value of 9 pulsars in a 10$\degr$ radius on average.
Our Monte Carlo simulations gave a probability of chance alignment smaller
than 0.01.

These calculations are back-of-the-envelope and cannot be the basis to claim
that Calvera and the radio ring are aligned, but they are a good indication
that it is very unlikely that their positional coincidence is by chance.

##  4 The nature of the radio ring

###  4.1 An H II region

H II regions have a radio spectral index of $>-0.2$ (Condon & Ransom, 2016),
which is already quite different from the best-fit spectral index value of
$-0.71\pm 0.09$ that we find here. H II regions also emit profusely in
H$\alpha$ (e.g., Haffner et al., 2009), which we do not detect. Moreover, the
area around Calvera was observed as part of the Wide-field Infrared Survey
Explorer (WISE, Wright et al., 2010) all-sky survey. We examined the four WISE
bands for emission coincident with the radio ring and detected no emission
above the image noise. Anderson et al. (2014) argue that the WISE survey
should be able to detect the mid-infrared emission from all Galactic H II
regions. Moreover, Makai et al. (2017) calculated the correlations between the
radio and IR flux densities of Galactic H II regions. Using the values (for
regions $>1$ pc) from their table 2, and a radio spectral index of $-0.1$, we
would expect $F_{8\leavevmode\nobreak\ \mu\mathrm{m}}\sim 40$ Jy,
$F_{12\leavevmode\nobreak\ \mu\mathrm{m}}\sim 80$ Jy and
$F_{22\leavevmode\nobreak\ \mu\mathrm{m}}\sim 90$ Jy if the ring were an H II
region. Even for a source with an area $>1,800$ arcmin2, these values are well
within the extended source sensitivities of their respective WISE bands.

In addition to the electromagnetic evidence against the H II interpretation,
the morphology also points against it. The circular shape of the source seems
to imply that it is ionised from within; however, a query of a 1 degree radius
around the coordinates of Calvera in the Catalog of Galactic OB Stars (Reed,
2003) yielded no results. We therefore discard the possibility that the ring
is an H II region.

###  4.2 A supernova remnant

Figure 3: 1 GHz surface brightness versus diameter for the Vukotić et al.
(2019) sample of Galactic SNRs with known distances (blue) and the Urošević et
al. (2005) sample of extragalactic SNRs. The purported Calvera SNR is plotted
in red using the measured surface brightness at 144 MHz and the best-fit
spectral index of 0.71, and the 54 pc diameter derived from the 3.3 kpc
distance estimate of Mereghetti et al. (2021).

The measured radio spectral index of $-0.71\pm 0.09$ already presents a good
indication that the radio emission is of a synchrotron origin, characteristic
of SNRs (Dubner & Giacani, 2015). In addition to this, there are two main
pieces of evidence that support a SNR interpretation for the ring around
Calvera: the low probability of a chance alignment between the pulsar and the
radio source (see section 3.3), and the presence of the extended X-ray source
1RXS J140818.1+792113. The fact that a possible association exists between all
three objects is what strengthens the SNR interpretation. Of course, there
still exists the possibility that the source is a SNR with no relation to the
Calvera pulsar other than the positional overlap. However, this would require
two stars exploding in the same line-of-sight, at a Galactic latitude where
SNRs are infrequent (the Green, 2019, catalogue lists 294 SNRs, only three of
which are at $b>10\degr$, and none are above $b=15\degr$). Our favoured
interpretation is that this radio ring is the remnant of the supernova
explosion that formed the Calvera pulsar. This would make the source SNR
G118.4+37.0, or Calvera’s SNR.

Perhaps the most convincing piece of evidence that the radio emission is from
a SNR is 1RXS J140818.1+792113: Zane et al. (2011) find extended soft X-ray
emission in XMM-Newton MOS and ROSAT All-Sky Survey (RASS) maps (indicated
with a dashed line in Figs. 1a and 1d). They perform a spectral analysis on
the data and find that the spectrum is well-fitted by a non-equilibrium
ionisation (NEI) model with an overabundance of oxygen; the former is a
typical signature of SNRs in general, and the latter of SNRs with a core-
collapse origin (Vink, 2012). The extended emission is also visible in the
public Chandra HRC-I observations used in Halpern & Gotthelf (2015). The
region of extended emission has angular separation of 13′ with Calvera, which
lead Zane et al. (2011) to suggest that the two are unrelated. However, both
are encompassed by the ring, and all three (the pulsar, the NEI plasma, and
the shell) could well be the remains of a single stellar explosion. The
additional possibility that only two of these components are related – the
plasma and the ring but not the pulsar, or any such combination – exists.
However, the presence of a SNR plasma inside the radio ring is very indicative
that the ring is a SNR, irrespective of whether the two are associated with
Calvera. A final point to be made regarding the NEI plasma is that the
relation between the H$\alpha$ smudge and the X-ray emitting material remains
unclear. Follow-up spectroscopy of the smudge could determine whether both
have similar compositions.

Mereghetti et al. (2021) estimated the distance to the Calvera pulsar to be
3.3 kpc (by assuming a neutron star radius of 13 km, and fitting the thermal
emission), and that the supernova explosion occurred in the Galactic halo. At
this distance, the SNR would have a diameter of 54 pc, and it would be at a
height over the Galactic plane of 2 kpc. This would make Calvera’s SNR, along
with the remnant of SN1006 (Raymond et al., 2017), one of the known remnants
with a largest height above the Galactic plane (with the important difference
that Calvera’s SNR would be the remnant of a core-collapse explosion). Getting
the (massive) progenitor star to such a large height above the Galactic plane
could require a high velocity runaway scenario: a star with a vertical
velocity of 100 km s-1 would take $\sim 30$ Myrs to reach it, which points to
a star in the low-mass end of core-collapse explosions
($8-10\leavevmode\nobreak\ \,{\rm M}_{\odot}$), or alternatively a star in
some kind of binary interaction that extends its lifetime (e.g., as discussed
in Zapartas et al., 2021). Hypervelocity runaway stars are rare (Renzo et al.,
2019) but certainly can be responsible for core-collapse supernovae in the
outskirts of galaxies (e.g., Zinn et al., 2011). Alternatively, the neutron
star need not be formed by a massive star at all, if its origin is the
accretion-induced or merger-induced collapse of a white dwarf (Grindlay &
Bailyn, 1988; Saio & Nomoto, 1985). This solves the problem of the location of
the progenitor, although the oxygen overabundance noted in the X-ray plasma by
Zane et al. (2011) favours a core-collapse origin.

For SNRs, there exists a relationship between their radio surface brightness
at 1 GHz and their diameters, known as the $\Sigma-D$ relation. Some authors
(e.g., Green, 2004) are sceptical of the soundness of this relation; however,
when we plot the diameter as derived from the Mereghetti et al. (2021)
distance against the radio surface brightness for the purported Calvera SNR
(see Fig. 3), the data point falls roughly in line with the samples of
Galactic and extragalactic SNRs found in Vukotić et al. (2019) and Urošević et
al. (2005), respectively. Using the values for the $\Sigma-D$ relation as
found in Vukotić et al. (2019), for a 1 GHz surface brightness of $1.74\times
10^{-23}$ W m2 Hz-1 sr-1 (using the best-fit spectral index of $-0.71$ and the
area between the inner and outer radii), we find a diameter of 62 pc, quite
comparable to the 54 pc diameter as derived from the Mereghetti et al. (2021)
distance.

An X-ray proper motion measurement by Halpern & Gotthelf (2015) places the
location of Calvera at the time of explosion $\sim 20$′ away from the center
of the ring (for an explosion that happened 8,000 years ago, see discussion
below; for an explosion corresponding to the characteristic age of the pulsar,
this distance is $>6$°). This could be interpreted as an indication that the
two are not related. However, the uncertainties in the Chandra observations
that produced this result (dominated by the short 2 ks HRC-I exposure obtained
in 2007), and the uncertainty in the age of the pulsar itself, are large.
Moreover, an inhomogeneous environment can easily produce an offset between
the SN explosion site and the SNR’s geometric center (Dohm-Palmer & Jones,
1996). This need not result in an asymmetric remnant: Williams et al. (2013)
showed that a uniform explosion into a density gradient produces a circular
remnant, but the center of the circle is different from the explosion site.
These are plausible reason for the observed misalignment; we cannot rule out
association simply from these proper motion measurements.

Another issue that draws attention to itself is that the X-ray plasma does not
overlap with the radio emission. This is uncommon, but not unheard of, in
SNRs: for instance, RCW 86 has a patch of X-ray emission with no radio
counterpart (Bamba et al., 2000). The presence of the optical smudge is also
puzzling: the only interpretation we can think of is that the SNR shock
encountered a small density discontinuity in its surrounding medium that
caused it to become radiative at that location.

Another source of tension comes from the pulsar characteristic age,
$\tau_{\mathrm{c}}=2.9\times 10^{5}$ yr (Halpern et al., 2013), which is well
past the time when we expect an SNR to still be visible. However,
characteristic ages often overestimate the true age of the pulsar; this is due
to the fact that the characteristic age neglects the initial period of the
pulsar. Since 59 ms is well within the distribution of birth periods inferred
from modeling the pulsar population (e.g. Faucher-Giguère & Kaspi, 2006)
Calvera’s characteristic age is only an upper limit to its true age. A case on
point are central compact objects (CCOs): thermal X-ray sources with no
accompanying radio or gamma-ray emission that are associated with known SNRs,
and which have characteristic ages $\sim 10^{8}$ yrs (De Luca, 2017).
Moreover, although in general SNRs live for $\sim 50,000$ years, we do not
know what is the spread in lifetimes –it could be intrinsically large, given
that a source’s lifetime is likely highly dependent on its environmental
conditions. An SNR in the halo, where the ambient densities are low, could
well live longer than 50,000 years. However, the Sedov-Taylor self-similar
solution (Sedov, 1959; Taylor, 1950), for an explosion energy of 1051 erg and
a SNR radius of 54 pc gives a SNR age much younger than that: $7,700$ years
for an ambient density of $4\times 10^{-4}$ cm-3 (we estimated the ambient
density in the halo around Calvera using the model of Miller & Bregman, 2013,
with $r=9.5$ kpc the distance between Calvera and the center of the Galaxy).
We do not know whether the conditions for self-similar evolution are present
on this remnant, though, so this age estimate, too, needs to be taken with a
grain of salt.

Finally, the ring shows no clear optical line emission. Abundant optical line
emission is a common feature in aged SNRs: when a SNR shock is slowed down to
approximately 200 km s-1 it begins to emit profusely in optical forbidden
lines, and cools quickly (Vink, 2012). However, a SNR’s cooling time and
length scales depend on the pre-shock density (Raymond, 1979), and in low
density environments the SNR enters the radiative phase at lower velocities
(corresponding to older ages). It is possible that the optical line emission
is too faint for our optical observations to recover, or that the remnant
still has not reached its radiative phase.

###  4.3 An Odd Radio Circle

Odd Radio Circles, or ORCs (Norris et al., 2021b), are a recently discovered
class of diffuse extragalactic objects, first detected in the Pilot Survey of
the Evolutionary Map of the Universe (Norris et al., 2021c). They are faint,
arcminute-sized circles of steep-spectrum radio emission at high Galactic
latitudes with no corresponding diffuse emission at optical, infrared, or
X-ray wavelengths (Norris et al., 2021a). As of this writing, there are only
seven published ORCs (Norris et al., 2021a; Koribalski et al., 2021; Filipović
et al., 2022; Omar, 2022), three of which have an elliptical galaxy at their
centre (referred to as their ‘host’ galaxy). It is unclear what phenomena
produce ORCs, and there is no evidence that the radio shells are associated
with their ‘host’ galaxies other than the positional coincidence of the two.
In fact, the most recently reported ORC (Filipović et al., 2022) was suggested
to be the remnant of an exploded runaway star that left its galaxy (the Large
Magellanic Cloud), whereas the remaining ones are thought to be distant
($z>0.1$) objects.

If further work disproves that the ring is a SNR, then it could be LOFAR’s
second contribution to the list of known ORCs (the first being Omar, 2022).
The Pan-STARRS PS1 images and catalog (Chambers et al., 2016b) shows no
prominent galaxy near the center of the ring as seen in other ORCs; out of the
18 PS1 objects located within $1\arcmin$ from the center of the ring, only a
single source (WISE J141113.52+792336.4) appears slightly extended
($1\arcsec\times 2\arcsec$) and has WISE colors consistent with those of
spiral galaxies. It is unclear whether this could be the ORC host galaxy if
the radio ring would be identified as an ORC. However, the known ORCs are
arcminute-sized objects, whereas the ring is $\sim 1$ degree in diameter. If
the ring is an extra-galactic object with comparable redshift to the other
ORCs, this would require that the mechanism that produces the ORCs is able to
create structures that span 1.5 orders of magnitude. A more likely scenario is
that this is an ORC of the type found by Filipović et al. (2022): a nearby
source, rather than a distant object. A nearer location would solve the issue
of its large scale but would take us back to the SNR interpretation. A final
point is that, morphologically, the ORCs appear to be quite different from the
ring discussed in this paper, as they tend to have more structure and be less
symmetric.

##  5 Discussion

We are fairly confident that this ring of radio emission seen in LOFAR, WENSS,
and the NVSS is the remnant of the explosion that formed Calvera. The main
pieces of evidence are the presence of an X-ray plasma with enhanced metals
and the rarity of a chance alignment with a pulsar at such high Galactic
latitudes. We discard the H II region interpretation given the lack of
infrared and H$\alpha$ emission; the ORC interpretation remains a possibility,
but the difference in size with the other known ORCs and the lack of a clear
host galaxy candidate weaken it. An additional possibility is that the ring is
a SNR unrelated to Calvera, although this would require two supernova
explosions occurring in a single line-of-sight a high Galactic latitude, which
is unlikely.

More data are necessary to determine whether the ring is an SNR or an ORC. In
particular, X-ray spectroscopy of the brighter filament at the North West
could unambiguously settle whether the ring is indeed a SNR. We do not know if
the ring shows no optical emission at all, or whether the optical emission is
too faint for our INT observations to detect, but an optical follow-up of the
ring will require deep exposures on dark conditions. Another set of follow-up
observations that could be more interesting is optical spectroscopy of the
H$\alpha$ smudge, to clarify whether metals are also present in its spectrum,
and whether it has a similar composition to the X-ray emitting plasma.
Finally, LOFAR LBA observations of the area at 58 MHz will better constrain
the radio spectral index, and polarisation studies at GHz frequencies can also
help settle whether the ring is a SNR.

There are some implications for the Calvera pulsar if the radio ring is a SNR
associated with it. Several authors (Zane et al., 2011; Halpern, 2011; Halpern
et al., 2013) have suggested that Calvera might be an ‘orphaned’ or aged CCO,
whose associated SNR either had faded away or was too faint to detect. The
discovery of a SNR around it would place Calvera in the CCO category. The age
of the pulsar would almost certainly be different from its characteristic age,
since the remnant seems to be on the younger side.

If confirmed to be an SNR, the ring around Calvera would be one of few SNRs
with large heights above the Galactic plane. This would make it an interesting
object for studies of SNR evolution in diffuse environments. Moreover, SNRs
shock and heat the material they expand into, and in this way make visible the
medium around them. SNRs at high Galactic latitudes are in a unique position
to probe the interstellar medium (ISM) of the Milky Way halo. The combined
analysis of deep X-ray and radio surveys, such as the one conducted by the
eROSITA telescope (Predehl et al., 2021) onboard the SRG observatory (Sunyaev
et al., 2021), and the suite of surveys carried out with MeerKAT and LOFAR
could discover a population of low surface brightness, high Galactic latitude
remnants (e.g., Churazov et al., 2021, 2022; Becker et al., 2021, or the
Calvera SNR in this work). These remnants could be longer-lived and fainter
than SNRs in the Milky Way disk, given the low ambient densities in the halo.

##  6 Summary

  1. 1.

We report the LOFAR 144 MHz discovery of a ring of diffuse radio emission
centered at $\alpha=14\mathrm{h}11\mathrm{m}12\aas@@fstack{s}6$,
$\delta=+79\degr 23\arcmin 15\arcsec$, with an outer radius
$R_{\mathrm{out}}=29\arcmin$, an inner radius $R_{\mathrm{inn}}=14\arcmin$,
and an integrated flux density of $1.08\pm 0.15$ Jy. The ring center is offset
4.87′ from the location of the Calvera pulsar. The ring shows low-significance
emission at 325 MHz and 1.4 GHz, and no emission in either of the four WISE
bands. Its spectral index is $-0.71\pm 0.09$.

  2. 2.

We did not detect radio emission at 144 MHz from the Calvera pulsar; we report
a 225 $\mu$Jy 3$\sigma$ upper limit for Calvera’s flux density at 144 MHz. We
also fail to detect pulsed radio emission at the known spin period of Calvera,
or as single pulses down to flux density and fluence limits of 0.25 mJy and
2.6 Jy ms at 149 MHz, respectively.

  3. 3.

We detect no radiative emission from the ring in H$\alpha$. We do detect a
smudge of extended H$\alpha$ emission internal to the ring.

  4. 4.

The non-detection of the ring in WISE, and its radio spectral index value
allow us to discard the possibility that the ring is an H II region.

  5. 5.

We estimate the probability of chance alignment between the pulsar and the
radio ring to be $<1$%.

  6. 6.

We favour an interpretation whereby the ring is the remnant of the supernova
explosion that formed Calvera, although more evidence is needed to support
this scenario.

  7. 7.

A final possibility is that the ring is an Odd Radio Circle, but this is
disfavoured due to the ring’s difference in size with other known ORCs.

  8. 8.

If the ring is an SNR, it is uniquely suited to probe the ISM in the Galactic
halo.

###### Acknowledgements.

LOFAR (van Haarlem et al., 2013) is the LOw Frequency ARray designed and
constructed by ASTRON. It has observing, data processing, and data storage
facilities in several countries, which are owned by various parties (each with
their own funding sources), and are collectively operated by the ILT
foundation under a joint scientific policy. The ILT resources have benefitted
from the following recent major funding sources: CNRS-INSU, Observatoire de
Paris and Université d’Orléans, France; BMBF, MIWF-NRW, MPG, Germany; Science
Foundation Ireland (SFI), Department of Business, Enterprise and Innovation
(DBEI), Ireland; NWO, The Netherlands; The Science and Technology Facilities
Council, UK; Ministry of Science and Higher Education, Poland; Istituto
Nazionale di Astrofisica (INAF), Italy. This research made use of the Dutch
national e-infrastructure with support of the SURF Cooperative (e-infra
180169) and the LOFAR e-infra group. The Jülich LOFAR Long Term Archive and
the German LOFAR network are both coordinated and operated by the Jülich
Supercomputing Centre (JSC), and computing resources on the supercomputer
JUWELS at JSC were provided by the Gauss Centre for Supercomputing e.V. (grant
CHTB00) through the John von Neumann Institute for Computing (NIC). This
research made use of the University of Hertfordshire high-performance
computing facility and the LOFAR-UK computing facility located at the
University of Hertfordshire and supported by STFC [ST/P000096/1], and of the
Italian LOFAR IT computing infrastructure supported and operated by INAF, and
by the Physics Department of Turin University (under an agreement with
Consorzio Interuniversitario per la Fisica Spaziale) at the C3S Supercomputing
Centre, Italy. The Isaac Newton Telesccope is operated on the island of La
Palma by the Isaac Newton Group of Telescopes in the Spanish Observatorio del
Roque de los Muchachos of the Instituto de Astrofísica de Canarias. The WFC
imaging was obtained as part of the Leiden education program. MA acknowledges
support from the VENI research programme with project number 202.143, which is
financed by the Netherlands Organisation for Scientific Research (NWO). AB
acknowledges support from the VIDI research programme with project number
639.042.729, which is financed by the Netherlands Organisation for Scientific
Research (NWO). RJvW acknowledges support from the ERC Starting Grant
ClusterWeb 804208. MJH acknowledges support from the UK Science and Technology
Facilities Council [ST/V000624/1].

## References

  * Anderson et al. (2014) Anderson, L. D., Bania, T. M., Balser, D. S., et al. 2014, ApJS, 212, 1
  * Bamba et al. (2000) Bamba, A., Koyama, K., & Tomida, H. 2000, PASJ, 52, 1157
  * Becker et al. (2021) Becker, W., Hurley-Walker, N., Weinberger, C., et al. 2021, A&A, 648, A30
  * Bertin & Arnouts (1996) Bertin, E. & Arnouts, S. 1996, A&AS, 117, 393
  * Bogdanov et al. (2019) Bogdanov, S., Ho, W. C. G., Enoto, T., et al. 2019, ApJ, 877, 69
  * Brentjens & de Bruyn (2005) Brentjens, M. A. & de Bruyn, A. G. 2005, A&A, 441, 1217
  * Chambers et al. (2016a) Chambers, K. C., Magnier, E. A., Metcalfe, N., et al. 2016a, arXiv e-prints, arXiv:1612.05560
  * Chambers et al. (2016b) Chambers, K. C., Magnier, E. A., Metcalfe, N., et al. 2016b, arXiv e-prints, arXiv:1612.05560
  * Churazov et al. (2021) Churazov, E. M., Khabibullin, I. I., Bykov, A. M., et al. 2021, MNRAS, 507, 971
  * Churazov et al. (2022) Churazov, E. M., Khabibullin, I. I., Bykov, A. M., et al. 2022, MNRAS, 513, L83
  * Cohen et al. (2007) Cohen, A. S., Lane, W. M., Cotton, W. D., et al. 2007, AJ, 134, 1245
  * Condon et al. (1998) Condon, J. J., Cotton, W. D., Greisen, E. W., et al. 1998, AJ, 115, 1693
  * Condon & Ransom (2016) Condon, J. J. & Ransom, S. M. 2016, Essential Radio Astronomy
  * Cordes & Lazio (2002) Cordes, J. M. & Lazio, T. J. W. 2002, ArXiv e-prints [arXiv:0207156]
  * De Luca (2017) De Luca, A. 2017, in Journal of Physics Conference Series, Vol. 932, Journal of Physics Conference Series, 012006
  * Dohm-Palmer & Jones (1996) Dohm-Palmer, R. C. & Jones, T. W. 1996, ApJ, 471, 279
  * Dubner & Giacani (2015) Dubner, G. & Giacani, E. 2015, A&A Rev., 23, 3
  * Erben et al. (2005) Erben, T., Schirmer, M., Dietrich, J. P., et al. 2005, Astronomische Nachrichten, 326, 432
  * Faucher-Giguère & Kaspi (2006) Faucher-Giguère, C.-A. & Kaspi, V. M. 2006, ApJ, 643, 332
  * Filipović et al. (2022) Filipović, M. D., Payne, J. L., Alsaberi, R. Z. E., et al. 2022, MNRAS[arXiv:2201.10026]
  * Gaensler & Johnston (1995) Gaensler, B. M. & Johnston, S. 1995, MNRAS, 277, 1243
  * Green (2004) Green, D. A. 2004, Bulletin of the Astronomical Society of India, 32, 335
  * Green (2019) Green, D. A. 2019, Journal of Astrophysics and Astronomy, 40, 36
  * Grindlay & Bailyn (1988) Grindlay, J. E. & Bailyn, C. D. 1988, Nature, 336, 48
  * Haberl (2007) Haberl, F. 2007, Ap&SS, 308, 181
  * Haffner et al. (2009) Haffner, L. M., Dettmar, R. J., Beckman, J. E., et al. 2009, Reviews of Modern Physics, 81, 969
  * Halpern (2011) Halpern, J. P. 2011, ApJ, 736, L3
  * Halpern et al. (2013) Halpern, J. P., Bogdanov, S., & Gotthelf, E. V. 2013, ApJ, 778, 120
  * Halpern & Gotthelf (2015) Halpern, J. P. & Gotthelf, E. V. 2015, ApJ, 812, 61
  * Hessels et al. (2007) Hessels, J. W. T., Stappers, B. W., Rutledge, R. E., Fox, D. B., & Shevchuk, A. H. 2007, A&A, 476, 331
  * Hotan et al. (2004) Hotan, A. W., van Straten, W., & Manchester, R. N. 2004, PASA, 21, 302
  * Intema et al. (2017) Intema, H. T., Jagannathan, P., Mooley, K. P., & Frail, D. A. 2017, A&A, 598, A78
  * Kondratiev et al. (2016) Kondratiev, V. I., Verbiest, J. P. W., Hessels, J. W. T., et al. 2016, A&A, 585, A128
  * Koribalski et al. (2021) Koribalski, B. S., Norris, R. P., Andernach, H., et al. 2021, MNRAS, 505, L11
  * Lane et al. (2014) Lane, W. M., Cotton, W. D., van Velzen, S., et al. 2014, MNRAS, 440, 327
  * Lorimer et al. (2006) Lorimer, D. R., Faulkner, A. J., Lyne, A. G., et al. 2006, MNRAS, 372, 777
  * Makai et al. (2017) Makai, Z., Anderson, L. D., Mascoop, J. L., & Johnstone, B. 2017, ApJ, 846, 64
  * Manchester et al. (2005) Manchester, R. N., Hobbs, G. B., Teoh, A., & Hobbs, M. 2005, AJ, 129, 1993
  * Mdzinarishvili & Melikidze (2004) Mdzinarishvili, T. G. & Melikidze, G. I. 2004, A&A, 425, 1009
  * Mereghetti et al. (2021) Mereghetti, S., Rigoselli, M., Taverna, R., et al. 2021, ApJ, 922, 253
  * Miller & Bregman (2013) Miller, M. J. & Bregman, J. N. 2013, ApJ, 770, 118
  * Norris et al. (2021a) Norris, R. P., Crawford, E., & Macgregor, P. 2021a, Galaxies, 9, 83
  * Norris et al. (2021b) Norris, R. P., Intema, H. T., Kapińska, A. D., et al. 2021b, PASA, 38, e003
  * Norris et al. (2021c) Norris, R. P., Marvil, J., Collier, J. D., et al. 2021c, PASA, 38, e046
  * Offringa et al. (2014) Offringa, A. R., McKinley, B., Hurley-Walker, N., et al. 2014, MNRAS, 444, 606
  * Offringa & Smirnov (2017) Offringa, A. R. & Smirnov, O. 2017, MNRAS, 471, 301
  * Omar (2022) Omar, A. 2022, Research Notes of the American Astronomical Society, 6, 100
  * Predehl et al. (2021) Predehl, P., Andritschke, R., Arefiev, V., et al. 2021, A&A, 647, A1
  * Ransom (2001) Ransom, S. M. 2001, PhD thesis, Harvard University
  * Raymond (1979) Raymond, J. C. 1979, ApJS, 39, 1
  * Raymond et al. (2017) Raymond, J. C., Winkler, P. F., Blair, W. P., & Laming, J. M. 2017, ApJ, 851, 12
  * Reed (2003) Reed, B. C. 2003, AJ, 125, 2531
  * Rengelink et al. (1997) Rengelink, R. B., Tang, Y., de Bruyn, A. G., et al. 1997, A&AS, 124, 259
  * Renzo et al. (2019) Renzo, M., Zapartas, E., de Mink, S. E., et al. 2019, A&A, 624, A66
  * Rutledge et al. (2008) Rutledge, R. E., Fox, D. B., & Shevchuk, A. H. 2008, ApJ, 672, 1137
  * Saio & Nomoto (1985) Saio, H. & Nomoto, K. 1985, A&A, 150, L21
  * Schirmer (2013) Schirmer, M. 2013, ApJS, 209, 21
  * Sedov (1959) Sedov, L. I. 1959, Similarity and Dimensional Methods in Mechanics
  * Shimwell et al. (2022) Shimwell, T. W., Hardcastle, M. J., Tasse, C., et al. 2022, A&A, 659, A1
  * Shimwell et al. (2017) Shimwell, T. W., Röttgering, H. J. A., Best, P. N., et al. 2017, A&A, 598, A104
  * Shimwell et al. (2019) Shimwell, T. W., Tasse, C., Hardcastle, M. J., et al. 2019, A&A, 622, A1
  * Sunyaev et al. (2021) Sunyaev, R., Arefiev, V., Babyshkin, V., et al. 2021, A&A, 656, A132
  * Tasse et al. (2021) Tasse, C., Shimwell, T., Hardcastle, M. J., et al. 2021, A&A, 648, A1
  * Taylor (1950) Taylor, G. 1950, Proceedings of the Royal Society of London Series A, 201, 159
  * Urošević et al. (2005) Urošević, D., Pannuti, T. G., Duric, N., & Theodorou, A. 2005, A&A, 435, 437
  * van Haarlem et al. (2013) van Haarlem, M. P., Wise, M. W., Gunst, A. W., et al. 2013, A&A, 556, A2
  * van Straten & Bailes (2011) van Straten, W. & Bailes, M. 2011, PASA, 28, 1
  * van Weeren et al. (2021) van Weeren, R. J., Shimwell, T. W., Botteon, A., et al. 2021, A&A, 651, A115
  * Vink (2012) Vink, J. 2012, A&A Rev., 20, 49
  * Voges et al. (1999) Voges, W., Aschenbach, B., Boller, T., et al. 1999, A&A, 349, 389
  * Vukotić et al. (2019) Vukotić, B., Ćiprijanović, A., Vučetić, M. M., Onić, D., & Urošević, D. 2019, Serbian Astronomical Journal, 199, 23
  * Williams et al. (2013) Williams, B. J., Borkowski, K. J., Ghavamian, P., et al. 2013, ApJ, 770, 129
  * Wright et al. (2010) Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868
  * Yao et al. (2017) Yao, J. M., Manchester, R. N., & Wang, N. 2017, ApJ, 835, 29
  * Yusifov & Küçük (2004) Yusifov, I. & Küçük, I. 2004, A&A, 422, 545
  * Zane et al. (2011) Zane, S., Haberl, F., Israel, G. L., et al. 2011, MNRAS, 410, 2428
  * Zapartas et al. (2021) Zapartas, E., de Mink, S. E., Justham, S., et al. 2021, A&A, 645, A6
  * Zinn et al. (2011) Zinn, P. C., Grunden, P., & Bomans, D. J. 2011, A&A, 536, A103


# Dynamical Phonons Following Electron Relaxation Stages in Photo-excited
Graphene

Nina Girotto  Centre for Advanced Laser Techniques, Institute of Physics,
10000 Zagreb, Croatia  Dino Novko<EMAIL_ADDRESS>Centre for Advanced
Laser Techniques, Institute of Physics, 10000 Zagreb, Croatia

###### Abstract

Ultrafast electron-phonon relaxation dynamics in graphene hides many distinct
phenomena, such as hot phonon generation, dynamical Kohn anomalies, and phonon
decoupling, yet still remains largely unexplored. Here, we unravel intricate
mechanisms governing the vibrational relaxation and phonon dressing in
graphene at a highly non-equilibrium state by means of first-principles
techniques. We calculate dynamical phonon spectral functions and momentum-
resolved linewidths for various stages of electron relaxation and find photo-
induced phonon hardening, overall increase of relaxation rate and
nonadiabaticity as well as phonon gain. Namely, the initial stage of photo-
excitation is found to be governed by strong phonon anomalies of finite-
momentum optical modes along with incoherent phonon production. Population
inversion state, on the other hand, allows production of coherent and
strongly-coupled phonon modes. Our research provides vital insights into the
electron-phonon coupling phenomena in graphene, and serves as a foundation for
exploring non-equilibrium phonon dressing in materials where ordered states
and phase transitions can be induced by photo-excitation.

###### keywords:

phonon dynamics, electron-phonon coupling, graphene, density functional theory

Through the ionic motion manipulation, photoexcitation as in pump-probe setup
is a powerful tool which paves the way for a highly effective designing and
customizing the desired functionalities of materials 1, 2. Namely, it can
induce novel phases 3, sometimes unreachable in equilibrium, opening, for
example, the possibility of light-induced superconductivity 4, 5, charge-
density-wave order 6, 7, ferroelectricity 8, and disorder-assisted structural
transition 9. Very often, these states of matter are characterized with a
strongly-coupled phonon mode and considerable electron-phonon coupling (EPC),
which are believed to be additionally renormalized in a photo-excited non-
equilibrium regime 10, 11. For instance, photo-induced softening of the
relevant phonon mode is quite common in ultrafast dynamics 12, 13,
nonetheless, photo-excitation can in some cases lead to phonon hardening and
consequently stabilize the structural phase 14, 15, 16, 17, 18.

Graphene exhibits extraordinary mechanical, transport and optoelectronic
properties 19, and is therefore an ideal platform to investigate the
fundamentals of ultrafast electron-lattice dynamics. Experimental techniques
such as time- and angle-resolved photoemission spectroscopy (tr-ARPES) 20, 21,
22, 23, 24, two-photon photoemission 25, 26, 27, and transient optical
spectroscopy 28, 29, 30 have reveled various aspects of carrier thermalization
in graphene, such as rapid electron-electron and electron-phonon recombination
of highly non-equilibrium distribution towards population inversion 21,
scattering of electrons with strongly-coupled optical phonons 22, 31, cooling
of hot carriers via acoustic phonons 22, as well as electron band
renormalization 30. On the other hand, non-equilibrium phonon dynamics in
graphene, phonon relaxation pathways and the corresponding photo-induced
phonon dressing are less investigated. Raman and coherent phonon spectroscopy
have demonstrated considerable phonon hardening of the $E_{2g}$ optical phonon
mode 14, 15, 32, which was attributed to the reduction of the nonadiabatic
electron-phonon interaction in non-equilibrium 14. Recent attosecond core-
level spectroscopy also uncovered ultrafast phonon stiffening of both zone-
center $E_{2g}$ and zone-edge $A_{1}^{\prime}$ phonon Kohn anomalies 33.

The theoretical studies of the aforesaid ultrafast phenomena are mostly based
on the two- and multi-temperature models 22, 34, 35, as well as on the time-
dependent Boltzmann equations 29, 36, and were proven to be valuable in
comprehending energy transfer between electrons and strongly-coupled hot
phonons. However, these methods do not account for transient phonon
renormalization and as such are not suitable for exploring all aspects of EPC
and phonon dynamics far from equilibrium, such as structural transitions and
soft phonon physics. The phonon renormalization in graphene was recently
inspected by means of real-time time-dependent density functional theory in
combination with molecular dynamics and real-space lattice distortions, which
allows for time-resolved self-consistent renormalization of phonons and EPC
strengths 11. However, since it relies on real-space distortions within the
supercell approach it is able to track the phonon dynamics of only few
selected modes. In addition, the study delivered somehow conflicting results,
i.e., instead of phonon hardening and electron-phonon decoupling as observed
in the experiments 14, 15, 32, 33, the phonon softening and enhanced EPC
strength were reported 11. This calls for further theoretical insights on non-
equilibrium phonon dynamics in graphene.

Here, we overcome these difficulties and investigate the effects of the photo-
excited population on phonon dynamics and EPC in graphene by means of
constrained density functional perturbation theory (cDFPT) 37, 38. Important
advantage of this approach is that it provides a full momentum-dependent
picture of phonon renormalization due to constrained photo-excited electron
distribution 39, 40, 37, 41, 42, 43. In addition, we combine cDFPT with
nonadiabatic phonon self-energy calculations in order to provide information
on phonon relaxation rates (linewidths) and nonadiabatic frequency
modifications, which are absent in the standard adiabatic cDFPT studies. We
discuss phonon renormalization for the usual stages of carrier relaxation in
graphene, namely, for strong non-equilibrium, population inversion, and hot
(Fermi-Dirac) carrier distribution. We observe remarkable modifications of the
well-known Kohn anomalies at the $\Gamma$ and K points of the Brillouin zone,
as well as appearance of the additional phonon anomalies away from the
Brillouin zone center induced by non-equilibrium population, renormalizing the
equilibrium dispersion for up to 6 meV. Light-induced increase of the overall
phonon linewidth and nonadiabatic effects are observed along with a striking
phonon gain, where the latter becomes coherent once graphene reaches the state
of photo-inversion. From the fermiology analysis we show that the EPC coupling
matrix elements are slightly reduced in non-equilibrium state, while the
observed features mostly stem from the modified scattering phase space under
transient conditions. With this work, we expand our microscopic understanding
of phonon relaxation dynamics of graphene in far non-equilibrium, which along
with the well-explored electron thermalization paths constitutes the full
dynamical picture of electron-phonon scatterings in graphene.

Figure 1: Schematic representation of different stages of electron relaxation
in photo-excited graphene. (a) Dirac cone with an equilibrium Fermi-Dirac
electron distribution. (b) At $t_{0}$, laser pulse excites electrons form the
[$\varepsilon_{1}$, $\varepsilon_{2}$] energy interval vertically upwards in
the conduction band, where they fill the states in the [$\varepsilon_{3}$,
$\varepsilon_{4}$] range. (c) Immediately after the pulse, electrons scatter
with other electrons and strongly-coupled phonons ($E_{2g}\simeq 200$ meV and
$A^{\prime}_{1}\simeq 160$ meV) until the population inversion is created. (d)
When the time scale of electron thermalization time ($\tau_{\rm th}$) is
reached, electrons follow a hot Fermi-Dirac distribution, which in our
calculations, amounts to 2200K.

Photo-excitation implies promoting a certain carrier density form the valence
to the conduction band separated by the energy of laser pulse $\hbar\Omega$,
which typically amounts to 1.5 eV. In this work, we track phonon properties
for various electron distributions inside the Dirac cone following directly
the pulse application. First, an equilibrium distribution [Fig. 1 (a)] is
disrupted with a short (fs) pulse causing an occurrence of empty states below
the Dirac point and filled states above it [Fig. 1 (b)]. Photo-generated
electrons and holes establish separate distributions and begin a rapid process
of thermalization and cooling through carrier-carrier and carrier-phonon
scatterings. The energy transfer to the strongly-coupled optical phonons
produces hot phonons, reported to exist on a femtosecond time-scale 28, 33.
Then, the photo-inverted state is established [Fig. 1 (c)] through the
competition between phonon-induced intraband scattering and Auger
recombination 44. The formation of the population inversion has already been
thoroughly explored with tr-ARPES 21, 45, 22, which reveals its relaxation
time of $\sim$100 fs. After its decay electrons follow a Fermi-Dirac
distribution at elevated temperatures [Fig. 1 (d)]. The whole process of
electron thermalization conceptually follows the one described for graphite in
Ref. 20. Subsequent hot-carrier cooling is governed by phonon-assisted
scatterings on the time scale of $1-10$ ps.

Figure 2: Adiabatic DFPT phonon dispersion in the presence of a photo-excited
(red line) and photo-inverted (yellow line) electron distribution in
comparison with an equilibrium result (grey dashed line). Two insets are a
zoom in region around high-symmetry points ($\Gamma$ and K), showing the
hardening of strongly-coupled optical modes ($E_{2g}$ and $A^{\prime}_{1}$)
along with a schematics of the corresponding atomic motions. Figure 3:
Dynamical phonon spectral functions at different stages of electron relaxation
in comparison with the adiabatic equilibrium DFPT result (grey dashed line),
i.e., for: (a) Equilibrium regime, (b) far non-equilibrium following the laser
excitation, (c) population inversion, (d) hot equilibrium distribution. Note
the strong renormalizations occurring near the two strongly-coupled optical
modes ($E_{2g}$ and $A_{1}^{\prime}$). Also, the negative linewidth
contribution to spectral function is shown in teal. The histograms in the
insets of (b)-(d) show the nonadiabatic corrections to the $E_{2g}$ and
$A^{\prime}_{1}$ modes ($\Delta^{\rm NA}=\omega^{\rm NA}-\omega^{\rm A}$). (e)
Phonon linewidth on $\Gamma$ to K path of LO/TO optical modes, due to EPC.
Color-coding is the same as in Fig. 2. (f-g) Intra- and inter-valley (i.e.,
$\mathrm{K}\rightarrow\mathrm{K}$ and
$\mathrm{K}\rightarrow\mathrm{K^{\prime}}$) electron transitions which
contribute to the phonon self-energy. The color-coded arrows reveal positive
(brown) and negative (teal) contributions to the linewidth. (h) Same as in (e)
but for the photo-inverted (yellow) and hot electron distributions (dark red).

We implemented the nonequlibrium distributions in the calculation of
electronic-structure properties and then calculated the renormalization of
phonons with the PHonon46, 47 and EPW 48, 49 codes in the adiabatic and
nonadiabatic approximations (see Supporting Information (SI) for more
details). The resulting adiabatic phonon dispersions are shown in Fig. 2. We
compare the phonon dispersion for the case of photo-excited electron
population as in Fig. 1(b) and photo-inverted distribution Fig. 1(c) with the
equilibrium adiabatic case. The largest effects of the nonequilibrium electron
distribution happen for the strongly-coupled optical modes $E_{2g}$ and
$A^{\prime}_{1}$ around $\Gamma$ and K points, respectively. We observe phonon
hardening for both phonons, and it ranges from $\simeq 1$ meV for the $E_{2g}$
mode and $4-6$ meV for the $A^{\prime}_{1}$ mode. Our results are in a good
agreement with Refs. 14, 15, 32 and especially with Ref. 33 where the
attosecond core-level spectroscopy demonstrated that the Raman-inactive
$A^{\prime}_{1}$ mode is the dominating channel for dissipation of electronic
coherence due to stronger coupling to electrons 50. Note also that in our
calculations the acoustic sum rule is not fully fulfilled because we inhibit
the long-wavelength dipole-allowed transitions 51, 52 by filling the states
above the Fermi energy.

In Fig. 3 we present the results of the non-equilibrium EPC calculations as
obtained with cDFPT. They contain the phonon spectral function
$B^{c}_{\nu}(\mathbf{q},\omega)$, which incorporates nonadiabatic
renormalization effects and transverse and longitudinal optical phonon
linewidth along the $\Gamma-\mathrm{K}$ path. With the term “adiabatic” we
refer to a calculation where the phonon energy $\omega$ is omitted in the
phonon self-energy calculations, and by “nonadiabatic” or “dynamic” where it
is included 53, 54 (see also SI). The first row represents four spectral
functions corresponding to the four distinct electron distributions from Fig.
1, together with the equilibrium adiabatic result for comparison. Equilibrium
graphene linewidth contributions for the $E_{2g}$ and $A^{\prime}_{1}$ modes
come from the vertical interband electron scattering inside the Dirac cone
($q\simeq\Gamma$) or between the two neighboring ones ($q\simeq\mathrm{K}$).
The features around $\Gamma$ and K points are symmetrical, but differ due to
the disparate EPC strengths.

The photo-excited electron distribution opens up new scattering possibilities,
which are schematically shown with arrows in Fig. 3(f). Besides the
significant modifications of the well-known nonadiabatic Kohn anomalies 55, 56
at the $\Gamma$ and K points coming from the non-equilibrium population, the
spectral function shows additional new anomalies further away from the
$\Gamma$ and K points. These photo-induced dynamical phonon anomalies come
from the electron transitions away from the Dirac point, at the photo-doped
regions. Compared to the equilibrium adiabatic dispersions, a 4 meV
renormalization is visible directly in the $\Gamma$ point for the $E_{2g}$
mode, and away from it, the highest optical branch is renormalized by 5 meV.
For the $A^{\prime}_{1}$ mode, we observe a 5 meV hardening and a 6 meV
modification at the intersection of the two optical branches. We note that
these sharp transient frequency modifications are quite large and are
comparable to the nonadiabatic frequency shifts of the $E_{2g}$ mode in the
highly-doped graphene 56. In the case of population inversion, the non-
equilibrium electron distribution is condensed in the vicinity of the Dirac
point, bringing the non-equilibrium spectral features closer to the $\Gamma$
and K points. We again observe phonon renormalization for the $E_{2g}$ and
$A^{\prime}_{1}$ modes of about 2 meV for both. Interestingly, for the strong
non-equilibrium and population inversion we observe additional phonon
hardening (softening) for $E_{2g}$ ($A_{1}^{\prime}$) when the nonadiabatic
effects are taken into account [insets of Figs. 3(b) and 3(c)]. In fact,
significant increase of the nonadiabatic correction is obtained for the strong
non-equilibrium, while it is reduced for the population inversion and almost
diminished for the hot equilibrium case. Note that this is contrary to the
conclusions drawn in Ref. 14, where the decrease of the nonadiabaticity is
suggested.

The corresponding contributions to the linewidth [Figs. 3(e) and 3(h)] show
that values at the $\Gamma$ and K points are unaltered, while slightly away
from these points it is significantly enhanced compared to its value at
equilibrium. For highly non-equilibrium case, additional notable phonon
broadening arise displaced from the high-symmetry points, at momenta where the
new dynamical anomalies appear. As stated, these linewidth features stem from
the electron transitions between the photo-excited block of filled states
above the Dirac point, and empty block below it.

Figure 4:  (a) Non-equilibrium phonon renormalization for different densities
of excited electrons. In the first row, we show the corresponding adiabatic
DFPT dispersions in the vicinity of $\Gamma$ and K points and observe phonon-
hardening increase with photo-excited electron density. The second row
contains the corresponding EPC induced linewidths. (b) Two cases of electron-
doped photo-excited graphene (i.e., for $E_{F}=0.2$ eV and $E_{F}=0.4$ eV). In
the first row, we show the corresponding dynamic spectral functions. The
second row shows the linewidths for the two doping regimes. We observe the
same features as in Fig. 3(e) with significant increase in the linewidth close
to the $\Gamma$ point. It derives from the larger electron density at the
Fermi level. Again, note the negative linewidth occurrence and dynamical
anomalies.

A crucial thing to notice is that in the photo-excited state, electrons can
scatter from the filled states at higher energies to the low-energy empty
states [see Fig. 3 (f-g), downwards pointing teal arrows], causing a negative
linewidth contribution. This phonon gain, happens in the immediate vicinity of
dynamical anomaly. For the population inversion, the phonon-gain contributions
are located directly at the $\Gamma$ and K high-symmetry points. In graphene,
acoustic phonon generation was experimentally achieved 57 and theoretically
explained 58. A conceptually similar phenomenon has recently been widely
explored in graphene, namely the photo-induced plasmon amplification with a
high potential for the development of novel optoelectronic devices 59, 60, 61,
62, 63, 64, 65, 66, 67, 68. Our observation of phonon gain is in agreement
with the observation of the negative plasmon linewidth and negative
conductivity in the photo-inverted state, and it simply means that hot phonons
are emitted in the non-equilibrium regime. In particular, the results show
that far non-equilibrium state supports the generation of incoherent hot
phonons with momenta slightly away from the high symmetry points (i.e., phonon
displacement pattern is shifted in phase between neighbouring unit cells),
while the population inversion supports coherent phonon generation of hot
$E_{2g}$ and $A_{1}^{\prime}$ (i.e., phonon displacement pattern is repeated
between neighbouring unit cells). This could explain, on the one hand, the
reduction of the phonon dephasing rate of $E_{2g}$ mode as reported in Ref.
14, and, on the other hand, the non-displacive mechanism for the generation of
both $E_{2g}$ and $A_{1}^{\prime}$ hot coherent phonons as obtained in
attosecond core-level spectroscopy 33.

As for the hot electron distribution, we calculated the number of carriers
located above the Dirac point in the state of photo-inversion and found the
temperature for which a Fermi-Dirac distribution produces the same number of
carriers in the conduction band. We show the spectral function for the hot
equilibrium electron distribution at $T=2200$ K. The $E_{2g}$ and
$A^{\prime}_{1}$ modes are slightly hardened, and here one can clearly see
also the edge of the electron-hole pair excitation continuum. The linewidth
resembles the one obtained for the photo-inverted population, only without the
negative contributions.

Further analysis includes changing the excited carrier density, which
experimentally corresponds to changing the laser fluence (Fig. 4). Effective
carrier density is denoted in Fig. 4(a) above the Dirac cones, and is
calculated as the summed density of photo-excited electrons and holes. As
expected, we observe larger phonon stiffening in the DFPT calculation for a
larger carrier density. Here we show only the adiabatic dispersions to focus
solely on the increased phonon stiffening in the $\Gamma$ and K points. The
$E_{2g}$ phonon hardening increases by 2 meV as the photo-excited density
increases by $6\times 10^{12}$ cm-1. Further, we observe larger phonon
linewidths deriving from the modified phase space with increasing carrier
density. Directly in the $\Gamma$ point the linewidth remains at its
equilibrium value, while slight differences in the position of the peaks away
from the high-symmetry points are visible. Furthermore, since in experiments
graphene is frequently placed on a substrate, we also provide results for
doped graphene, specifically for Fermi levels of $E_{F}=200$ and 400 meV. We
compare the photo-doped spectral function with the adiabatic DFPT equilibrium
calculation for the corresponding Fermi energy (dashed black lines). Again, we
notice larger phonon hardening around the K point (5 meV for both dopings)
then around the $\Gamma$ point (4 meV for $E_{F}=200$ meV and 1 meV for
$E_{F}=400$ meV). We observe dynamical anomalies at the same positions as they
occur in the pristine photo-doped case. We notice how with increased doping,
the dynamic phonon dispersion softens and the effects of photo-induced phonon
hardening are less pronounced. Largest softening is observed around the
dynamical anomalies. In general, for doped graphene, the intrinsic linewidths
of the two highest optical modes are larger then for pristine graphene. When
doped graphene is photo-excited, linewidth behaves in the same fashion as in
the pristine photo-doped case, with the strong dynamical anomalies occurring
in the vicinity of the high-symmetry points.

Figure 5: The analysis of static phonon self-energies, revealing the
electronic processes behind phonon anomalies in the case of photo-excited
(first row) and photo-inverted (second-row) electron distribution. (a)
Relative change in the $\mathbf{k}$-summed EPC matrix elements
$\Delta^{rel}|g_{\nu}(\mathbf{k})|^{2}$ for the two highest optical modes and
along the $\Gamma-\mathrm{K}$ q-path. (b)
$\Delta^{rel}|g_{E_{2g}}(\mathbf{k})|^{2}$ resolved in $\mathbf{k}$ space.
(c),(d) $\mathbf{k}$-resolved static susceptibility $\chi^{c}_{0}(\mathbf{k})$
for the $E_{2g}$ and $A^{\prime}_{1}$ modes, respectively. The positive
(negative) contributions to both $\Delta^{rel}|g_{E_{2g}}(\mathbf{k})|^{2}$
and $\chi^{c}_{0}(\mathbf{k})$ are shown with red (blue).

Finally, we present the analysis of the adiabatic phonon self-energy as
calculated in cDFPT $\pi^{c}_{\nu}(\mathbf{q})$ 69, 70 (see also SI).
Averaging out the electronic degrees of freedom deriving from the EPC matrix
elements, leads to
$|g_{\nu}^{nm,c}(\mathbf{k},\mathbf{q})|\rightarrow|g_{\nu}^{c}(\mathbf{q})|$
and the self-energy expression can be written as
$\pi^{c}_{\nu}(\mathbf{q})=|g^{c}_{\nu}(\mathbf{q})|^{2}\chi^{c}_{0}(\mathbf{q})$,
where $\chi^{c}_{0}(\mathbf{q})$ denotes the bare susceptibility function.

In this way, we can separate the non-equilibrium effects that come from the
modifications in screened EPC matrix elements $|g^{c}_{\nu}(\mathbf{q})|$ and
the photo-induced changes in the available phase space via
$\chi^{c}_{0}(\mathbf{k})$. In Fig. 5 we show the analysis for the conduction
band, but the results for the valence band are the same due to the electron-
hole symmetry. In the first two columns we show the relative difference
between the EPC matrix elements in the photo-excited and photo-inverted
graphene with respect to the pristine equilibrium case for the $E_{2g}$ and
$A_{1}^{\prime}$ modes
$\Delta^{rel}|g_{\nu}(\mathbf{k})|^{2}=(|g^{c}_{\nu}(\mathbf{k})|^{2}-|g^{\rm
eq}_{\nu}(\mathbf{k})|^{2})/{|g^{\rm eq}_{\nu}(\mathbf{k})|^{2}}$. In the
first column, the relative differences are calculated after the electron-
phonon matrix elements were summed throughout the whole Brillouin zone for a
chosen $\mathbf{q}$ point on the $\Gamma-\mathrm{K}$ path and for each one of
the two highest optical modes. The largest value of the relative change is
only $\pm$10% and it appears for those wavevectors $\mathbf{q}$ for which the
electronic transitions are forbidden by the specific electronic structure of
graphene and for which the EPC strength is weak at equilibrium (i.e., away
from the $\Gamma$ and K points). More dramatic changes occur if
$\Delta^{rel}|g_{E_{2g}}(\mathbf{k})|^{2}$ is resolved in the kx-ky plane
[column (b) in Fig. 5]. Here we explicitly show the result for
$\mathbf{q}=\Gamma$ for which we found the largest values of
$\Delta^{rel}|g_{E_{2g}}(\mathbf{k})|^{2}$ (see SI for results in additional q
points). For both the photo-excited and photo-inverted electron distribution
case, our calculations show that $\Delta^{rel}|g_{E_{2g}}(\mathbf{k})|^{2}$
reaches values of $\pm$ 40% in certain regions of kx-ky space for the E2g
mode. We observe a symmetrical pattern of positive and negative contributions
to $\Delta^{rel}|g_{E_{2g}}(\mathbf{k})|^{2}$. Due to phase space
restrictions, only a small region around the K points is picked out when doing
a self-energy calculation, or multiplying $|g^{c}_{\nu}(\mathbf{q})|^{2}$ and
$\chi^{c}_{0}(\mathbf{q})$. In other words, as shown in columns (c) and (d) of
Fig. 5, $\chi^{c}_{0}(\mathbf{k})$ is finite around the Dirac points also in a
symmetric pattern. In this way, when calculating the phonon self-energy and
summing over the whole Brillouin zone, the net effect of the large relative
changes $\Delta^{rel}|g_{E_{2g}}(\mathbf{k})|^{2}$ cancels out as equal
amounts of positive and negative values are picked out by the
$\chi^{c}_{0}(\mathbf{k})$ factor (see Sec. S2 in SI for more detailed
discussion).

Therefore, the effects of photo-induced changes of the electron-phonon matrix
elements are small in the end, and it turns out that the phonon hardening and
decrease of the dephasing rate observed in Ref. 14, along with other non-
equilibrium phenomena presented here, come dominantly from the photo-induced
changes in the carrier scattering phase space $\chi^{c}_{0}(\mathbf{k})$.
Therefore, optically-induced phonon blueshifts are not necessarily a sign of
the suppressed coupling and it could come from the pure electronic origin, as
was shown for instance in the case of photo-excited TiSe2 18. This resolves
the debate on whether the EPC strength is suppressed 14 or enhanced 11 in the
non-equilibrium graphene.

We elaborate this claim by thoroughly inspecting the color-coded contributions
to the bare electron susceptibility $\chi^{c}_{0}(\mathbf{k})$. The features
for $\mathbf{q}=\mathrm{K}$ and $\mathbf{q}=\Gamma$ are in essence the same,
so the discussion is applicable to both. In the equilibrium case where the
Fermi surface is almost a point, the only contribution comes from the Dirac
point. Since the electrons can only scatter from the filled states below the
Fermi energy to the empty states above it, the final results for
$\chi^{0}(\mathbf{k})$ is negative. Focusing now on the photo-excited case
(Fig. 5, first row), we first notice the mentioned equilibrium contribution
(red dots) positioned directly in the Dirac points. Photo-excited electrons
fill the states visible here as an additional triangle around each Dirac
point. Each triangle consists of positive and negative susceptibility
contributions. Electrons from the positive $\chi^{c}_{0}(\mathbf{k})$ portion
are responsible for the negative linewidth contribution, as they contribute by
scattering from the higher energy filled states to the lower energy empty
ones. Due to finite temperature, this is in principle also possible within an
equilibrium distribution, but those contributions are then suppressed by the
much larger negative $\chi_{0}(\mathbf{k})$ contribution. The electrons from
the negative $\chi^{c}_{0}(\mathbf{k})$ section make standard transitions to
higher energy empty states. In the EPC calculations, we used a dynamical
susceptibility term [see Eq. (S1) in SI], which, together with varying the
wavevector $\mathbf{q}$, leads to the competition between these two
contributions and, hence, to the net negative linewidth regions visible in
Fig. 3. The obvious difference between susceptibility contribution shape in
the $\Gamma$ and K points, is the result of trigonal warping 71, which is
reversed in the neighboring K points. Setting $\mathbf{q}=\mathrm{K}$ leads to
the superposition of two relatively rotated triangles, making the
$A^{\prime}_{1}$ susceptibility look circular. $\chi^{c}_{0}(\mathbf{k})$ for
the photo-inverted distribution, consists of two circular contributions, as
the filled/empty states are in the energy range of a linear electron
distribution. Color-coding again suggests that the electrons closer to the
Dirac point rather scatter to the empty states below the Fermi level, while
the higher energy electrons do the opposite and contribute positively to the
phonon linewidth. In this case, varying the wavevector $\mathbf{q}$ in the
dynamical susceptibility calculation reduces the phase space for the electrons
in the negative section, as a number of vertical transitions is restricted. It
confines the phonon gain to be placed directly in the high-symmetry points.

Note in the end that the present analysis on phonon dynamics based in cDFPT is
not only restricted to the non-equilibrium induced by laser pulses, but it
could be utilized to study phonon renormalization in any out-of-equilibrium
conditions. For instance, impact of the static electric fields and the
corresponding non-equilibrium electrons on phonons in current-carrying
materials is still not well understood despite its importance for
comprehending various transport properties 72, 73, 74.

Understanding the mechanisms behind the out-of-equilibrium EPC provides
valuable insights into the fundamental physics of graphene and helps unravel
the complex interplay between charge carriers and lattice vibrations. We
investigated the coupling of the high-energy optical phonon modes with the
photo-excited electron distribution by means of cDFPT. We observed hardening
of the well-known Kohn anomalies at the center and edge of the Brillouin zone.
The latter comes from the modified phase space for electron transitions which
leads to different screening effects, while the effective EPC coupling
strenghts are only slightly changed. We obtained complex nonadiabatic EPC
features that emerge in both the dispersion and linewidth, and mostly
originate from the new scattering channels opened in non-equilibrium. For
instance, sharp dynamical phonon anomalies away from the high-symmetry points
and the overall increase of the phonon scattering rate have been observed.
Also, we showed incoherent phonon gain at finite wavevectors irrespective of
the doping level or the concentration of the photo-excited carriers, while
coherent phonon generation is expected in the state of population inversion
and is within the scope of typical experiments. We believe our work offers
crucial information on the nature of EPC in graphene, justifies the known non-
equilibrium features and sheds new light on the understanding of the
underlying ultrafast vibrational relaxation mechanisms.

Useful discussions with Jan Berges and Samuel Poncé are gratefully
acknowledged. We acknowledge financial support from the Croatian Science
Foundation (Grant no. UIP-2019-04-6869) and from the European Regional
Development Fund for the “Center of Excellence for Advanced Materials and
Sensing Devices” (Grant No. KK.01.1.1.01.0001).

More information on computational details and detailed analysis of the
modifications in the electron-phonon matrix elements and scattering phase
space due to non-equilibrium distribution

## References

  * Först et al. 2011 Först, M.; Manzoni, C.; Kaiser, S.; Tomioka, Y.; Tokura, Y.; Merlin, R.; Cavalleri, A. Nonlinear phononics: A new ultrafast route to lattice control. _Nat. Phys._ 2011, _7_
  * de la Torre et al. 2021 de la Torre, A.; Kennes, D. M.; Claassen, M.; Gerber, S.; McIver, J. W.; Sentef, M. A. Colloquium: Nonthermal pathways to ultrafast control in quantum materials. _Rev. Mod. Phys._ 2021, _93_ , 041002
  * Basov et al. 2017 Basov, D.; Averitt, R.; Hsieh, D. Towards properties on demand in quantum materials. _Nat. Mater._ 2017, _16_ , 1077–1088
  * Fausti et al. 2011 Fausti, D.; Tobey, R. I.; Dean, N.; Kaiser, S.; Dienst, A.; Hoffmann, M. C.; Pyon, S.; Takayama, T.; Takagi, H.; Cavalleri, A. Light-Induced Superconductivity in a Stripe-Ordered Cuprate. _Science_ 2011, _331_ , 189
  * Mitrano et al. 2016 Mitrano, M.; Cantaluppi, A.; Nicoletti, D.; Kaiser, S.; Perucchi, A.; Lupi, S.; Pietro, P. D.; Pontiroli, D.; Riccò, M.; Clark, S. R.; Jaksch, D.; Cavalleri, A. Possible light-induced superconductivity in K3C60 at high temperature. _Nature_ 2016, _530_ , 461
  * Stojchevska et al. 2014 Stojchevska, L.; Vaskivskyi, I.; Mertelj, T.; Kusar, P.; Svetin, D.; Brazovskii, S.; Mihailovic, D. Ultrafast Switching to a Stable Hidden Quantum State in an Electronic Crystal. _Science_ 2014, _344_ , 177
  * Kogar et al. 2020 Kogar, A. et al. Light-induced charge density wave in LaTe3. _Nat. Phys._ 2020, _16_ , 159
  * Nova et al. 2019 Nova, T. F.; Disa, A. S.; Fechner, M.; Cavalleri, A. Metastable ferroelectricity in optically strained SrTiO3. _Science_ 2019, _364_ , 1075
  * Wall et al. 2018 Wall, S. E.; Yang, S.; Vidas, L.; Chollet, M.; Glownia, J. M.; Kozina, M. E.; Katayama, T.; Henighan, T.; Jiang, M. P.; Miller, T. A.; Reis, D. A.; Boatner, L. A.; Delaire, O.; Trigo, M. Ultrafast disordering of vanadium dimers in photoexcited VO2. _Science_ 2018, _362_ , 572 – 576
  * Pomarico et al. 2017 Pomarico, E.; Mitrano, M.; Bromberger, H.; Sentef, M. A.; Al-Temimy, A.; Coletti, C.; Stöhr, A.; Link, S.; Starke, U.; Cacho, C.; Chapman, R.; Springate, E.; Cavalleri, A.; Gierz, I. Enhanced electron-phonon coupling in graphene with periodically distorted lattice. _Phys. Rev. B_ 2017, _95_ , 024304
  * Hu et al. 2022 Hu, S.-Q.; Zhao, H.; Lian, C.; Liu, X.-B.; Guan, M.-X.; Meng, S. Tracking photocarrier-enhanced electron-phonon coupling in nonequilibrium. _npj Quantum Mater._ 2022, _7_ , 14
  * Hase et al. 2002 Hase, M.; Kitajima, M.; Nakashima, S.-i.; Mizoguchi, K. Dynamics of Coherent Anharmonic Phonons in Bismuth Using High Density Photoexcitation. _Phys. Rev. Lett._ 2002, _88_ , 067401
  * Thiemann et al. 2022 Thiemann, F.; Sciaini, G.; Kassen, A.; Hagemann, U.; Meyer zu Heringdorf, F.; Horn-von Hoegen, M. Ultrafast transport-mediated homogenization of photoexcited electrons governs the softening of the ${A}_{1g}$ phonon in bismuth. _Phys. Rev. B_ 2022, _106_ , 014315
  * Ishioka et al. 2008 Ishioka, K.; Hase, M.; Kitajima, M.; Wirtz, L.; Rubio, A.; Petek, H. Ultrafast electron-phonon decoupling in graphite. _Phys. Rev. B_ 2008, _77_ , 121402
  * Yan et al. 2009 Yan, H.; Song, D.; Mak, K. F.; Chatzakis, I.; Maultzsch, J.; Heinz, T. F. Time-resolved Raman spectroscopy of optical phonons in graphite: Phonon anharmonic coupling and anomalous stiffening. _Phys. Rev. B_ 2009, _80_ , 121403
  * Krapivin et al. 2022 Krapivin, V. et al. Ultrafast Suppression of the Ferroelectric Instability in ${\mathrm{KTaO}}_{3}$. _Phys. Rev. Lett._ 2022, _129_ , 127601
  * Jiang et al. 2016 Jiang, M. P. et al. The origin of incipient ferroelectricity in lead telluride. _Nat. Commun._ 2016, _7_
  * Otto et al. 2021 Otto, M. R.; Pöhls, J.-H.; de Cotret, L. P. R.; Stern, M. J.; Sutton, M.; Siwick, B. J. Mechanisms of electron-phonon coupling unraveled in momentum and time: The case of soft phonons in TiSe2. _Sci. Adv._ 2021, _7_ , eabf2810
  * Basov et al. 2014 Basov, D. N.; Fogler, M. M.; Lanzara, A.; Wang, F.; Zhang, Y. Colloquium: Graphene spectroscopy. _Rev. Mod. Phys._ 2014, _86_ , 959–994
  * Rohde et al. 2018 Rohde, G.; Stange, A.; Müller, A.; Behrendt, M.; Oloff, L.-P.; Hanff, K.; Albert, T. J.; Hein, P.; Rossnagel, K.; Bauer, M. Ultrafast Formation of a Fermi-Dirac Distributed Electron Gas. _Phys. Rev. Lett._ 2018, _121_ , 256401
  * Gierz et al. 2013 Gierz, I.; Petersen, J. C.; Mitrano, M.; Cacho, C.; Turcu, I. C. E.; Springate, E.; Stöhr, A.; Köhler, A.; Starke, U.; Cavalleri, A. Snapshots of non-equilibrium Dirac carrier distributions in graphene. _Nat. Mater._ 2013, _12_ , 1119–1124
  * Johannsen et al. 2013 Johannsen, J. C.; Ulstrup, S.; Cilento, F.; Crepaldi, A.; Zacchigna, M.; Cacho, C.; Turcu, I. C. E.; Springate, E.; Fromm, F.; Raidel, C.; Seyller, T.; Parmigiani, F.; Grioni, M.; Hofmann, P. Direct View of Hot Carrier Dynamics in Graphene. _Phys. Rev. Lett._ 2013, _111_ , 027403
  * Na et al. 2019 Na, M. X.; Mills, A. K.; Boschini, F.; Michiardi, M.; Nosarzewski, B.; Day, R. P.; Razzoli, E.; Sheyerman, A.; Schneider, M.; Levy, G.; Zhdanovich, S.; Devereaux, T. P.; Kemper, A. F.; Jones, D. J.; Damascelli, A. Direct determination of mode-projected electron-phonon coupling in the time domain. _Science_ 2019, _366_ , 1231
  * Düvel et al. 2022 Düvel, M.; Merboldt, M.; Bange, J. P.; Strauch, H.; Stellbrink, M.; Pierz, K.; Schumacher, H. W.; Momeni, D.; Steil, D.; Jansen, G. S. M.; Steil, S.; Novko, D.; Mathias, S.; Reutzel, M. Far-from-Equilibrium Electron–Phonon Interactions in Optically Excited Graphene. _Nano Lett._ 2022, _22_ , 4897
  * Xu et al. 1996 Xu, S.; Cao, J.; Miller, C. C.; Mantell, D. A.; Miller, R. J. D.; Gao, Y. Energy Dependence of Electron Lifetime in Graphite Observed with Femtosecond Photoemission Spectroscopy. _Phys. Rev. Lett._ 1996, _76_ , 483–486
  * Moos et al. 2001 Moos, G.; Gahl, C.; Fasel, R.; Wolf, M.; Hertel, T. Anisotropy of Quasiparticle Lifetimes and the Role of Disorder in Graphite from Ultrafast Time-Resolved Photoemission Spectroscopy. _Phys. Rev. Lett._ 2001, _87_ , 267402
  * Tan et al. 2017 Tan, S.; Argondizzo, A.; Wang, C.; Cui, X.; Petek, H. Ultrafast Multiphoton Thermionic Photoemission from Graphite. _Phys. Rev. X_ 2017, _7_ , 011004
  * Kampfrath et al. 2005 Kampfrath, T.; Perfetti, L.; Schapper, F.; Frischkorn, C.; Wolf, M. Strongly Coupled Optical Phonons in the Ultrafast Dynamics of the Electronic Energy and Current Relaxation in Graphite. _Phys. Rev. Lett._ 2005, _95_ , 187403
  * Breusing et al. 2009 Breusing, M.; Ropers, C.; Elsaesser, T. Ultrafast Carrier Dynamics in Graphite. _Phys. Rev. Lett._ 2009, _102_ , 086809
  * Pagliara et al. 2011 Pagliara, S.; Galimberti, G.; Mor, S.; Montagnese, M.; Ferrini, G.; Grandi, M. S.; Galinetto, P.; Parmigiani, F. Photoinduced $\pi$-$\pi^{\ast}$ Band Gap Renormalization in Graphite. _J. Am. Chem. Soc._ 2011, _133_ , 6318
  * Kampfrath et al. 2005 Kampfrath, T.; Perfetti, L.; Schapper, F.; Frischkorn, C.; Wolf, M. Strongly Coupled Optical Phonons in the Ultrafast Dynamics of the Electronic Energy and Current Relaxation in Graphite. _Phys. Rev. Lett._ 2005, _95_ , 187403
  * Ferrante et al. 2018 Ferrante, C.; Virga, A.; Benfatto, L.; Martinati, M.; Fazio, D.; Sassi, U.; Fasolato, C.; Ott, A.; Postorino, P.; Yoon, D.; Cerullo, G.; Mauri, F.; Ferrari, A.; Scopigno, T. Raman spectroscopy of graphene under ultrafast laser excitation. _Nat. Commun._ 2018, _9_
  * Sidiropoulos et al. 2021 Sidiropoulos, T. P. H. et al. Probing the Energy Conversion Pathways between Light, Carriers, and Lattice in Real Time with Attosecond Core-Level Spectroscopy. _Phys. Rev. X_ 2021, _11_ , 041060
  * Caruso et al. 2020 Caruso, F.; Novko, D.; Draxl, C. Photoemission signatures of nonequilibrium carrier dynamics from first principles. _Phys. Rev. B_ 2020, _101_ , 035128
  * Caruso and Novko 2022 Caruso, F.; Novko, D. Ultrafast dynamics of electrons and phonons: from the two-temperature model to the time-dependent Boltzmann equation. _Adv. Phys.: X_ 2022, _7_ , 2095925
  * Tong and Bernardi 2021 Tong, X.; Bernardi, M. Toward precise simulations of the coupled ultrafast dynamics of electrons and atomic vibrations in materials. _Phys. Rev. Res._ 2021, _3_ , 023072
  * Murray et al. 2007 Murray, E. D.; Fahy, S.; Prendergast, D.; Ogitsu, T.; Fritz, D. M.; Reis, D. A. Phonon dispersion relations and softening in photoexcited bismuth from first principles. _Phys. Rev. B_ 2007, _75_ , 184301
  * Nomura and Arita 2015 Nomura, Y.; Arita, R. Ab initio downfolding for electron-phonon-coupled systems: Constrained density-functional perturbation theory. _Phys. Rev. B_ 2015, _92_ , 245108
  * Liu et al. 2022 Liu, K.; Mao, S.; Zhang, S.; Zhou, J. Photoinduced Rippling of Two-Dimensional Hexagonal Nitride Monolayers. _Nano Lett._ 2022, _22_ , 9006–9012
  * Tangney and Fahy 2002 Tangney, P.; Fahy, S. Density-functional theory approach to ultrafast laser excitation of semiconductors: Application to the ${A}_{1}$ phonon in tellurium. _Phys. Rev. B_ 2002, _65_ , 054302
  * Murray et al. 2005 Murray, E. D.; Fritz, D. M.; Wahlstrand, J. K.; Fahy, S.; Reis, D. A. Effect of lattice anharmonicity on high-amplitude phonon dynamics in photoexcited bismuth. _Phys. Rev. B_ 2005, _72_ , 060301
  * Paillard et al. 2019 Paillard, C.; Torun, E.; Wirtz, L.; Íñiguez, J.; Bellaiche, L. Photoinduced Phase Transitions in Ferroelectrics. _Phys. Rev. Lett._ 2019, _123_ , 087601
  * Marini and Calandra 2021 Marini, G.; Calandra, M. Lattice dynamics of photoexcited insulators from constrained density-functional perturbation theory. _Phys. Rev. B_ 2021, _104_ , 144103
  * Winzer et al. 2013 Winzer, T.; Malić, E.; Knorr, A. Microscopic mechanism for transient population inversion and optical gain in graphene. _Phys. Rev. B_ 2013, _87_ , 165413
  * Gierz et al. 2015 Gierz, I.; Mitrano, M.; Petersen, J. C.; Cacho, C.; Turcu, I. C. E.; Springate, E.; Stöhr, A.; Köhler, A.; Starke, U.; Cavalleri, A. Population inversion in monolayer and bilayer graphene. _J. Phys.: Condens. Matter_ 2015, _27_ , 164204
  * Giannozzi et al. 2017 Giannozzi, P. et al. Advanced capabilities for materials modelling with Quantum ESPRESSO. _J. Phys. Condens. Matter_ 2017, _29_ , 465901
  * Baroni et al. 2001 Baroni, S.; de Gironcoli, S.; Dal Corso, A.; Giannozzi, P. Phonons and related crystal properties from density-functional perturbation theory. _Rev. Mod. Phys._ 2001, _73_ , 515
  * Poncé et al. 2016 Poncé, S.; Margine, E.; Verdi, C.; Giustino, F. EPW: Electron–phonon coupling, transport and superconducting properties using maximally localized Wannier functions. _Comput. Phys. Commun._ 2016, _209_ , 116
  * Lee et al. 2023 Lee, H. et al. Electron–phonon physics from first principles using the EPW code. _npj Comput. Mater._ 2023, _9_ , 156
  * Krylow et al. 2020 Krylow, S.; Hernandez, F. V.; Bauerhenne, B.; Garcia, M. E. Ultrafast structural relaxation dynamics of laser-excited graphene: Ab initio molecular dynamics simulations including electron-phonon interactions. _Phys. Rev. B_ 2020, _101_ , 205428
  * van Loon et al. 2021 van Loon, E. G. C. P.; Berges, J.; Wehling, T. O. Downfolding approaches to electron-ion coupling: Constrained density-functional perturbation theory for molecules. _Phys. Rev. B_ 2021, _103_ , 205103
  * Berges et al. 2022 Berges, J.; Girotto, N.; Wehling, T.; Marzari, N.; Poncé, S. Phonon self-energy corrections: To screen, or not to screen. _arXiv_ 2022, _2212.11806_
  * Engelsberg and Schrieffer 1963 Engelsberg, S.; Schrieffer, J. R. Coupled Electron-Phonon System. _Phys. Rev._ 1963, _131_ , 993
  * Giustino 2017 Giustino, F. Electron-phonon interactions from first principles. _Rev. Mod. Phys._ 2017, _89_ , 015003
  * Piscanec et al. 2004 Piscanec, S.; Lazzeri, M.; Mauri, F.; Ferrari, A. C.; Robertson, J. Kohn Anomalies and Electron-Phonon Interactions in Graphite. _Phys. Rev. Lett._ 2004, _93_ , 185503
  * Lazzeri and Mauri 2006 Lazzeri, M.; Mauri, F. Nonadiabatic Kohn Anomaly in a Doped Graphene Monolayer. _Phys. Rev. Lett._ 2006, _97_ , 266407
  * Andersen et al. 2019 Andersen, T.; Dwyer, B.; Sanchez-Yamagishi, J.; Rodriguez-Nieva, J.; Agarwal, K.; Watanabe, K.; Taniguchi, T.; Demler, E.; Kim, P.; Park, H.; Lukin, M. Electron-phonon instability in graphene revealed by global and local noise probes. _Science_ 2019, _364_ , 154–157
  * Zhao et al. 2013 Zhao, C.; Xu, W.; Peeters, F. M. Cerenkov emission of terahertz acoustic-phonons from graphene. _Appl. Phys. Lett._ 2013, _102_ , 222101
  * Sabbaghi et al. 2015 Sabbaghi, M.; Lee, H.-W.; Stauber, T.; Kim, K. S. Drift-induced modifications to the dynamical polarization of graphene. _Phys. Rev. B_ 2015, _92_ , 195429
  * Morgado and Silveirinha 2022 Morgado, T. A.; Silveirinha, M. G. Directional dependence of the plasmonic gain and nonreciprocity in drift-current biased graphene. _Nanophotonics_ 2022, _11_ , 4929
  * Ryzhii et al. 2007 Ryzhii, V.; Ryzhii, M.; Otsuji, T. Negative dynamic conductivity of graphene with optical pumping. _J. Appl. Phys._ 2007, _101_ , 083114 – 083114
  * Mikhailov et al. 2016 Mikhailov, S. A.; Savostianova, N. A.; Moskalenko, A. S. Negative dynamic conductivity of a current-driven array of graphene nanoribbons. _Phys. Rev. B_ 2016, _94_ , 035439
  * Petersen et al. 2017 Petersen, R.; Pedersen, T. G.; Javier García de Abajo, F. Nonlocal plasmonic response of doped and optically pumped graphene, ${\mathrm{MoS}}_{2}$, and black phosphorus. _Phys. Rev. B_ 2017, _96_ , 205430
  * Morgado and Silveirinha 2017 Morgado, T. A.; Silveirinha, M. G. Negative Landau Damping in Bilayer Graphene. _Phys. Rev. Lett._ 2017, _119_ , 133901
  * Boubanga-Tombet et al. 2020 Boubanga-Tombet, S.; Knap, W.; Yadav, D.; Satou, A.; But, D. B.; Popov, V. V.; Gorbenko, I. V.; Kachorovskii, V.; Otsuji, T. Room-Temperature Amplification of Terahertz Radiation by Grating-Gate Graphene Structures. _Phys. Rev. X_ 2020, _10_ , 031004
  * Park et al. 2022 Park, S.; Sammon, M.; Mele, E.; Low, T. Plasmonic gain in current biased tilted Dirac nodes. _Nat. Commun._ 2022, _13_ , 7667
  * Low et al. 2018 Low, T.; Chen, P.-Y.; Basov, D. N. Superluminal plasmons with resonant gain in population inverted bilayer graphene. _Phys. Rev. B_ 2018, _98_ , 041403
  * Duppen et al. 2016 Duppen, B. V.; Tomadin, A.; Grigorenko, A. N.; Polini, M. Current-induced birefringent absorption and non-reciprocal plasmons in graphene. _2D Mater._ 2016, _3_ , 015011
  * Berges et al. 2020 Berges, J.; van Loon, E. G. C. P.; Schobert, A.; Rösner, M.; Wehling, T. O. Ab initio phonon self-energies and fluctuation diagnostics of phonon anomalies: Lattice instabilities from Dirac pseudospin physics in transition metal dichalcogenides. _Phys. Rev. B_ 2020, _101_ , 155107
  * Novko 2020 Novko, D. Broken adiabaticity induced by Lifshitz transition in MoS2 and WS2 single layers. _Commun. Phys._ 2020, _3_ , 1
  * Dresselhaus 1974 Dresselhaus, G. Graphite Landau levels in the presence of trigonal warping. _Phys. Rev. B_ 1974, _10_ , 3602–3609
  * Sabbaghi et al. 2022 Sabbaghi, M.; Stauber, T.; Lee, H.-W.; Gomez-Diaz, J. S.; Hanson, G. W. In-plane optical phonon modes of current-carrying graphene. _Phys. Rev. B_ 2022, _105_ , 235405
  * Mao et al. 2022 Mao, W.-H.; Shang, M.-Y.; Lü, J.-T. Hot and cold phonons in electrically biased graphene. _Phys. Rev. B_ 2022, _106_ , 125406
  * Mao et al. 2023 Mao, W.-H.; Shang, M.-Y.; Lü, J.-T. Current-driven collective dynamics of non-Hermitian edge vibrations in armchair graphene nanoribbons. _Phys. Rev. B_ 2023, _107_ , 085419