# Ionized Gas Extended Over 40 kpc in an Odd Radio Circle Host Galaxy

Alison L. Coil∗1 Serena Perrotta1 David S. N. Rupke2 Cassandra Lochhaas3
Christy A. Tremonti4 Aleks Diamond-Stanic5 Drummond Fielding6 James Geach7
Ryan C. Hickox8 John Moustakas9 Gregory H. Rudnick10 Paul Sell11 Kelly E.
Whalen8

###### Abstract

A new class of extragalactic astronomical sources discovered in 2021, named
Odd Radio Circles (ORCs)[1], are large rings of faint, diffuse radio continuum
emission spanning $\sim$1 arcminute on the sky. Galaxies at the centers of
several ORCs have photometric redshifts of $z\sim 0.3-0.6$, implying physical
scales of several 100 kiloparsecs in diameter for the radio emission, the
origin of which is unknown. Here we report spectroscopic data on an ORC
including strong [O ii] emission tracing ionized gas in the central galaxy of
ORC4 at $z=0.4512$. The physical extent of the [O ii] emission is $\sim$40 kpc
in diameter, larger than expected for a typical early-type galaxy[2] but an
order of magnitude smaller than the large-scale radio continuum emission. We
detect a $\sim$200 km s-1 velocity gradient across the [O ii] nebula, as well
as a high velocity dispersion of $\sim$180 km s-1. The [O ii] equivalent width
(EW, $\sim$50 Å) is extremely high for a quiescent galaxy. The morphology,
kinematics, and strength of the [O ii] emission are consistent with the infall
of shock ionized gas near the galaxy, following a larger-scale, outward moving
shock driven by a galactic wind. Both the extended optical and radio emission,
while observed on very different scales, may therefore result from the same
dramatic event.

Center for Astrophysics and Space Sciences, University of California, 9500
Gilman Dr., La Jolla, CA 92093<EMAIL_ADDRESS>
Department of Physics, Rhodes College, Memphis, TN 38112

Space Telescope Science Institute, 3700 San Martin Dr., Baltimore, MD 21218

Department of Astronomy, University of Wisconsin-Madison, Madison, WI 53706

Department of Physics and Astronomy, Bates College, Lewison, ME 04240

Center for Computational Astrophysics, Flatiron Institute, 162 5th Ave, New
York, NY 10010

Centre for Astrophysics Research, University of Hertfordshire, Hatfield,
Hertfordshire AL10 9AB

Department of Physics and Astronomy, Dartmouth College, Hanover, NH 03755

Department of Physics and Astronomy, Siena College, Loudonville, NY 12211

Department of Physics and Astronomy, University of Kansas, Lawrence, KS 66045

Department of Astronomy, University of Florida, Gainesville, FL, 32611

Three ORCs were first discovered[1] in the Evolutionary Map of the Universe
Pilot Survey[3] at a frequency of $\sim$1 GHz using the Australian Square
Kilometer Array Pathfinder (ASKAP) telescope[4, 5]. The ASKAP array has high
angular resolution ($\sim$13′′), is sensitive to low surface brightness
emission, and the pilot survey covers 270 deg2 of the sky, such that it is
able to detect rare, faint objects not previously detected. A fourth ORC
(ORC4) was discovered in archival data taken with the Giant MeterWave Radio
Telescope (GMRT)[6] at 325 MHz, while additional ORCs were discovered in later
ASKAP[7] and MeerKAT data[8]. ORCs exhibit large, limb-brightened rings of
radio continuum emission with lower surface brightness emission in the
interior. Follow-up MeerKAT data on ORC1 reveals that the radio emission is
most likely aged synchrotron, and polarization data show magnetic field lines
that are consistent with an expanding shell[9]. The origin of the radio
emission is unknown, with scenarios proposed including Galactic supernovae
remnants, double-lobed radio galaxies, ring galaxies, and gravitationally-
lensed Einstein rings[1]. However, the most likely scenarios involve a shock
from an outflowing galactic wind or a blast wave driven by merging
supermassive black holes[7].

To search for an optical counterpart to ORCs, we observed ORC4 with the Keck
Cosmic Web Imager (KCWI[10]) and detected strong [O ii] and weak Mg ii and [Ne
iii] line emission in a spatially-integrated spectrum (Extended Data Fig. 1).
The [O ii] luminosity is $7.1\times 10^{41}\ erg\ s^{-1}$, near the break of
the [O ii] luminosity function at $z\sim 0.5$[11]. The restframe [O ii] EW is
50 Å in the spatially-integrated spectrum, which is an order of magnitude
higher than what is typically found in red, early-type galaxies[12]. Long-slit
optical spectroscopy of ORC4 that includes wavelengths redward of the KCWI
data reveals low-ionization nuclear emission-line region (LINER)-like line
ratios (D. S. N. Rupke et al., manuscript in preparation).

A spectral energy distribution (SED) fit to the optical and infrared emission
in the central galaxy of ORC4 (Extended Data Fig. 2) reveals that the galaxy
is massive (log $M_{*}/M_{\odot}=11.27\pm 0.06$) and has an old stellar
population, with a stellar age of $t_{age}=6.0\ (+/-1.6)$ Gyr. While the SED
fit has only a 3% contribution from an AGN, the radio continuum emission in
this galaxy[1] is likely due to an AGN, as there is no indication of on-going
star formation. The observed radio flux of 1.43 mJy at 325 MHz corresponds to
a luminosity of 1.24$\times 10^{24}\ W\ Hz^{-1}$ at a restframe frequency of
200 MHz, given the redshift and radio spectral index[1]; this is 1.5 orders of
magnitude fainter than the break in the radio AGN luminosity function[13].

The [O ii] emission is bright enough to create spatially-resolved maps to
probe the morphology, extent, and resolved kinematics of the ionized gas. Fig.
1 shows a map of the [O ii] surface brightness, the peak of which is
$2.3\times 10^{-16}\ erg\ s^{-1}\ cm^{-2}\ arcsec^{-2}$. The azimuthally-
averaged surface brightness profile (Fig. 2) reveals a strong central
concentration of the [O ii] emission, in contrast to the larger-scale limb-
brightened radio emission. The radial extent is $\sim$20 kpc ($\sim$40 kpc in
diameter), and the half light (or “effective”) radius is 3.4 kpc, deconvolved
with the seeing. The [O ii] radial extent is an order of magnitude smaller
than the radio continuum emission in ORC4, which extends $\sim$200 kpc in
radius.

The [O ii] restframe EW spatial map (Fig. 1, panel d) does not show strong
asymmetry. The restframe EW at the center of the nebula is $>$70 Å(Fig. 2),
and at a radius of 15 kpc (shown here at $R/R_{e}=3$, relative to the galaxy
half-light radius) it is $\sim$25 Å. For comparison, [O ii] EW profiles are
shown for massive, nearby early-type galaxies, where integral field data
reveals that 38% of galaxies with $>10^{11.5}$ M⊙ have [O ii] emission[2]. The
median [O ii] EW of these comparison galaxies is 4.8 Å, and the maximum EW is
28 Ang. The [O ii] EW in ORC4 is therefore an order of magnitude higher than
is typically seen in massive, early-type galaxies. The maximal radial extent
of the [O ii] in these comparison galaxies varies from 0.6-18.2 kpc, with a
median maximal extent of 2.3 kpc. The [O ii] emission in ORC4, detected to a
radius of $\sim$20 kpc, is therefore unusually spatially extended. It is
extended not only in the sense that it is not from the nucleus of the galaxy,
but the maximal extent is an order of magnitude higher than in comparable
galaxies.

The kinematics of the [O ii] emission (Fig. 3) reveal an asymmetric velocity
gradient across the nebula, from $+$170 km s-1at the northern edge to $-50$ km
s-1at the southern edge. A north-south position-velocity diagram (Fig. 3,
panel c) highlights the asymmetry in the velocity gradient. The velocity
dispersion varies from a minimum of $\sim$60 km s-1 at the edges of the nebula
to $\sim$180 km s-1 in the center and up to $\sim 250$ km s-1 at the northern
and western edges.

There are multiple possibilities for the origin of the extended [O ii]
emission observed in ORC4. Given the properties of the ionized gas and host
galaxy, it is unlikely that it is associated with the interstellar medium
(ISM) of the galaxy or an AGN or AGN-driven outflow or due to a cooling flow
(see Methods section “Alternative scenarios for ionized gas origin” for
details.) Another possibility is that it is associated with the same energetic
event that created the large-scale radio emission which identified the source
as an ORC. The LINER-like line ratios of ORC4 point to the ionized gas arising
from shocks[14]. One of the possible scenarios for the creation of ORCs is a
shock from a starburst wind, which would create multiple shocks on different
scales including both a forward moving shock on large scales and a reverse
shock, sometimes called a wind shock, on smaller scales, with a contact
discontinuity in between[15, 16]. In a model of galactic winds created by
starburst activity[17], the forward shock can exist on scales of $>$100 kpc,
while the reverse shock exists on scales of $\sim$10-50 kpc. Once the
starburst episode that initially drove the outflow has shut off, the gas
heated by the reverse shock can then travel backwards towards the central
galaxy[18], moving to smaller scales. The shocked wind expands to fill the
under-pressurized inner region that has been cleared by the outward-moving
forward shock and may produce a turbulent, energetic medium full of additional
shocks by interacting with gas in and around the galaxy. Wind re-heated in
this way by the reverse shock may then radiatively cool giving rise to ionized
gas emission on the scale of the galaxy. Meanwhile the forward shock, which
has continued to propogate outwards, creates synchrotron emission on large
scales. This scenario differs than that proposed by Norris et al. (2022),
where the large-scale radio emission is associated with the reverse shock; we
propose that is associated with the forward shock. In theoretical models the
large-scale forward shock can propogate to scales of hundreds of kiloparsecs
on timescales of several Gyr[18]. Recent cosmological simulations focusing on
galaxy mergers as the cause of ORCs[19] show that additional shocks are seen
on smaller scales as well.

This scenario is consistent with observations of ORC4 in multiple ways. The
best fit SED model indicates that the galaxy had a burst of star formation
$\sim$1 Gyr ago and is no longer forming stars. The scale, morphology, and
surface brightness distribution of the [O ii] emission are consistent with
expectations for the shocked wind gas falling back towards the galaxy, which
could also lead to an asymmetric velocity gradient and high velocity
dispersion due to turbulence. The strength of the [O ii] emission is also
consistent with being due to shocks, which are known to lead to an enhancement
of [O ii][20], and/or may reflect mixing between hot and cool gas phases,
which can lead to additional gas cooling[21, 22].

To further test a starburst-driven wind model for ORCs, we ran a suite of
simulations that launch a wind from an isolated galaxy and followed the
evolution of the forward and reverse shocks (Fig. 4). Our goal was to
reproduce the forward shock radius ($\sim$200 kpc) and Mach number ($\sim$1-2)
required to produce the observed synchrotron emission, the scale of the
reverse shock and subsequent infalling shocked wind ($\sim$20 kpc), the
elevated velocity dispersion of the shocked wind gas ($\sim$150 km s-1), and
the age of the wind bubble ($\sim$1 Gyr). We found that these parameters could
be reproduced with a wind that is launched with a velocity of 450 km s-1 and a
mass outflow rate of 200 $M_{\odot}$yr-1, which blows for 200 Myr. (See
Methods section “Outflowing galactic wind simulation” for details.) The high
mass outflow rate needed is potentially consistent with the SED fit results
for the central galaxy in ORC4, which require a high fraction of the total
stellar mass (61%) produced in a starburst which presumably had a high SFR,
which could have created an outflow with a high mass outflow rate. The
features of the wind model presented here are generic features of an energy
and mass injection, without regard to the source of the injection. In order
for the wind bubble to survive at a Mach number of $\sim$1-2 to a scale of
$\sim$200 kpc, the density of the CGM must be lower than is typically assumed
for a galaxy of stellar mass $10^{11}$ $M_{\odot}$, though the density is
still feasible (simulations find the CGM density of a Milky Way-mass halo
spans six orders of magnitude[23]). As ORCs appear to be a rare phenomena,
they may arise only in situations where the galaxy progenitor and CGM
parameters are favorable. This simulation shows that a single event such as a
strong starburst can launch a wind that simultaneously produces both the
large-scale radio emission and the smaller-scale [O ii] emission scale, when
observed $\sim$1 Gyr after the event.

We have been studying outflowing galactic winds in massive ($\sim 10^{11}$ M⊙)
starburst galaxies at $0.4<z<0.8$[24]; most are late-stage major mergers[25].
These starburst galaxies are driving extremely fast multi-phase gas outflows
with velocities of $\sim$2000 km s-1 and mass outflow rates of $\gtrsim$200
$M_{\odot}$yr-1[26, 27]. Their space densities ($1\times 10^{-6}Mpc^{-3}$) are
similar to that of post-starburst galaxies and ultra-luminous infrared
galaxies (ULIRGs)[28]. KCWI data reveals [O ii] nebulae in these galaxies of
sizes 30-100 kpc[29] (S. Perrotta et al., in preparation).

Though their space density is potentially lower (estimated at $\sim 2\times
10^{-8}Mpc^{-3}$)[9], ORCs may be a later stage of a similar phenomenon. In
order for ORCs to be the result of a starburst wind shock the starburst must
have been extreme[9], with a high mass outflow rate, and the density of the
CGM must be low. The massive starburst galaxies we have been studying at
$z\sim 0.5$ have such extreme outflow properties[30]. In the model of an
ancient, starburst-driven wind, the ancestors of ORCs must resemble these
compact starbursts in some way. These ORC wind relics may then be a direct
tracer of powerful outflows from earlier energetic starburst events.

## References

  * [1] Norris, R. P. _et al._ Unexpected circular radio objects at high Galactic latitude.  _PASA_ 38, e003 (2021).
  * [2] Pandya, V. _et al._ The MASSIVE Survey. VI. The Spatial Distribution and Kinematics of Warm Ionized Gas in the Most Massive Local Early-type Galaxies.  _ApJ_ 837, 40 (2017).
  * [3] Norris, R. P. _et al._ The Evolutionary Map of the Universe pilot survey.  _PASA_ 38, e046 (2021).
  * [4] Johnston, S. _et al._ Science with the Australian Square Kilometre Array Pathfinder.  _PASA_ 24, 174 (2007).
  * [5] McConnell, D. _et al._ The Australian Square Kilometre Array Pathfinder: Performance of the Boolardy Engineering Test Array.  _PASA_ 33, e042 (2016).
  * [6] Ananthakrishnan, S.  The Giant Meterwave Radio Telescope / GMRT.  _Journal of Astrophysics and Astronomy Supplement_ 16, 427 (1995).
  * [7] Koribalski, B. S. _et al._ Discovery of a new extragalactic circular radio source with ASKAP: ORC J0102-2450.  _MNRAS_ 505, L11–L15 (2021).
  * [8] Koribalski, B. S. _et al._ MeerKAT discovery of a double radio relic and odd radio circle.  _arXiv e-prints_ arXiv:2304.11784 (2023).  2304.11784.
  * [9] Norris, R. P. _et al._ MeerKAT uncovers the physics of an odd radio circle.  _MNRAS_ 513, 1300–1316 (2022).
  * [10] Morrissey, P. _et al._ The Keck Cosmic Web Imager Integral Field Spectrograph.  _ApJ_ 864, 93 (2018).
  * [11] Zhu, G., Moustakas, J. & Blanton, M. R.  The [O II] $\lambda$3727 Luminosity Function at z ~1\.  _ApJ_ 701, 86–93 (2009).
  * [12] Yan, R. _et al._ On the Origin of [O II] Emission in Red-Sequence and Poststarburst Galaxies.  _ApJ_ 648, 281–298 (2006).
  * [13] Franzen, T. M. O. _et al._ The GLEAM 200-MHz local radio luminosity function for AGN and star-forming galaxies.  _PASA_ 38, e041 (2021).
  * [14] Yan, R. & Blanton, M. R.  The Nature of LINER-like Emission in Red Galaxies.  _ApJ_ 747, 61 (2012).
  * [15] Faucher-Giguère, C.-A. & Quataert, E.  The physics of galactic winds driven by active galactic nuclei.  _MNRAS_ 425, 605–622 (2012).
  * [16] King, A. & Pounds, K.  Powerful Outflows and Feedback from Active Galactic Nuclei.  _ARA &A_ 53, 115–154 (2015).
  * [17] Thompson, T. A., Quataert, E., Zhang, D. & Weinberg, D. H.  An origin for multiphase gas in galactic winds and haloes.  _MNRAS_ 455, 1830–1844 (2016).
  * [18] Lochhaas, C., Thompson, T. A., Quataert, E. & Weinberg, D. H.  Fast winds drive slow shells: a model for the circumgalactic medium as galactic wind-driven bubbles.  _MNRAS_ 481, 1873–1896 (2018).
  * [19] Dolag, K., Böss, L. M., Koribalski, B. S., Steinwandel, U. P. & Valentini, M.  Insights on the Origin of Odd Radio Circles from Cosmological Simulations.  _ApJ_ 945, 74 (2023).
  * [20] Allen, M. G., Groves, B. A., Dopita, M. A., Sutherland, R. S. & Kewley, L. J.  The MAPPINGS III Library of Fast Radiative Shock Models.  _ApJS_ 178, 20–55 (2008).
  * [21] Gronke, M. & Oh, S. P.  How cold gas continuously entrains mass and momentum from a hot wind.  _MNRAS_ 492, 1970–1990 (2020).
  * [22] Fielding, D. B. & Bryan, G. L.  The Structure of Multiphase Galactic Winds.  _ApJ_ 924, 82 (2022).
  * [23] Simons, R. C. _et al._ Figuring Out Gas & Galaxies in Enzo (FOGGIE). IV. The Stochasticity of Ram Pressure Stripping in Galactic Halos.  _ApJ_ 905, 167 (2020).  2004.14394.
  * [24] Tremonti, C. A., Moustakas, J. & Diamond-Stanic, A. M.  The Discovery of 1000 km $s^{-1}$ Outflows in Massive Poststarburst Galaxies at z=0.6.  _ApJ_ 663, L77–L80 (2007).
  * [25] Sell, P. H. _et al._ Massive compact galaxies with high-velocity outflows: morphological analysis and constraints on AGN activity.  _MNRAS_ 441, 3417–3443 (2014).
  * [26] Rupke, D. S. N. _et al._ The Ionization and Dynamics of the Makani Galactic Wind.  _ApJ_ 947, 33 (2023).
  * [27] Perrotta, S. _et al._ Kinematics, Structure, and Mass Outflow Rates of Extreme Starburst Galactic Outflows.  _ApJ_ 949, 9 (2023).  2303.07448.
  * [28] Whalen, K. E. _et al._ The Space Density of Intermediate-redshift, Extremely Compact, Massive Starburst Galaxies.  _AJ_ 164, 222 (2022).
  * [29] Rupke, D. S. N. _et al._ A 100-kiloparsec wind feeding the circumgalactic medium of a massive compact galaxy.  _Nature_ 574, 643–646 (2019).
  * [30] Diamond-Stanic, A. M. _et al._ Compact Starburst Galaxies with Fast Outflows: Central Escape Velocities and Stellar Mass Surface Densities from Multiband Hubble Space Telescope Imaging.  _ApJ_ 912, 11 (2021).
  * [31] Dey, A. _et al._ Overview of the DESI Legacy Imaging Surveys.  _AJ_ 157, 168 (2019).

Figure 1: Optical stellar continuum and [O ii] emission in ORC4. (a) DESI
Legacy Survey DR9[31] optical $grz$ image of ORC4, spanning 40′′ on each side.
The KCWI field of view (FOV) is shown with a yellow box, and the dotted white
box traces the central 10′′ $\times$ 10′′, shown in the remaining panels. (b)
Stellar continuum surface brightness image from KCWI, integrated over the
observed wavelength range 4200-5300 Å, with axes labeled in kpc from the
center of the galaxy. White contours show levels of 2.5, 5, 10, and 15
$\times\ 10^{-16}\ erg\ s^{-1}\ cm^{-2}\ arcsec^{-2}$. (c) [O ii] surface
brightness image, showing spaxels with S/N$\geq$3\. White contours show
surface brightness levels of 0.3, 0.8, 1.4, and 1.9 $\times\ 10^{-16}\ erg\
s^{-1}\ cm^{-2}\ arcsec^{-2}$. (d) [O ii] restframe equivalent width (EW) map.

Figure 2: Radial surface brightness and equivalength width profile of the [O
ii] emission. (a) The black line shows the [O ii] radial surface brightness
profile, with 1$\sigma$ errors; the nebula has a radial extent of $\sim$20
kpc. The green line shows the radial profile of the stellar continuum,
normalized to the peak of the [O ii] profile, and the orange line shows the
extent of the ground-based seeing, also normalized to the [O ii]. The red line
is the best fit slope to the [O ii] radial profile (beyond 5 kpc, outside the
seeing disk) of $-2.8$. (b) The black line shows the restframe equivalent
width (EW) profile of the [O ii] emission in ORC4 as a function of radius from
the center of the galaxy, relative to the stellar half light radius
($R_{e}=4.8$ kpc), with 1$\sigma$ errors. For comparison, purple dot-dashed
lines show the [O ii] EW radial profile fits for nearby massive, early-type
galaxies[2].

Figure 3: Kinematics of the [O ii] emission.  (a) Velocity map of the [O ii]
emission showing the central velocity ($V_{50}$) of the ionized gas relative
to the systemic redshift of the source, $z=0.4512$. Contours mark where the
central velocity is -50, -25, 0, 25, 75 km s-1. A clear velocity gradient is
seen from north to south across the nebula. (b) The velocity dispersion
$\sigma$ of the [O ii] emission, showing a higher dispersion in the center,
north, and west regions of the nebula. Insets show spectra at various
locations around the nebula, highlighting areas with high and low velocity
dispersion. Blue and cyan lines in the insets show Gaussian fits to the
individual [O ii] emission lines. (c) A position-velocity diagram through the
center of the [O ii] nebula along the north-south direction, averaging across
the central 5 columns, with 1$\sigma$ errors. (d) The velocity dispersion
through the center of the nebula along the north-south direction, averaging
across the central 5 columns, with 1$\sigma$ errors.

Figure 4: Simulation of starburst-driven wind.  (a) Two-dimensional maps of
gas temperature (upper left half of each panel) and gas radial velocity
relative to the galaxy (lower right half of each panel), shown such that red
colors indicate material falling onto the galaxy and blue colors indicate
material outflowing from the galaxy, at three times since the onset of the
starburst. The wind blows for 200 Myr, after which wind material that has been
heated by the reverse shock falls back towards the galaxy, while the forward
shock continues to expand into the halo. (b) Radius of the large-scale forward
shock and the smaller-scale inner and outer boundaries of the shocked wind as
a function of time since the onset of the starburst. (c) Spherical two-
dimensional map (mirrored from the one quadrant simulated) of the gas
temperature at the time of observation, 740 Myr after the onset of the
starburst and 540 Myr since the star formation and wind ceased. The large
yellow ring corresponds to where radio synchrotron emission is expected to
arise, while the smaller-scale black region surrounding the central galaxy
contains cool gas that could be observed as [O ii] emission.

##  0.1 KCWI observations and data reduction.

We observed ORC4 (RA=15:55:24.63 Dec=$+$27:26:34.3) with KCWI on the Keck II
telescope on March 27, 2022. The weather was clear and the seeing was
$\sim$0.9′′; we observed at a position angle of 0 degrees and an airmass of
1.2. We used the blue low dispersion (BL) grating and the medium slicer on
KCWI, which provides a spectral resolution of R$=$1800, a spaxel size of
0.29′′ $\times$ 0.69′′, and a FOV of 16′′ $\times$ 20′′ per pointing. We
configured the instrument to a central wavelength of 4700 Å and used a
detector binning of 2 $\times$ 2, which provides wavelength coverage of 3700
$-$ 5700 Å. We observed a total of three KCWI pointings, each spanning 16′′
$\times$ 20′′, integrating for 20 minutes per pointing. One pointing was
centered on ORC4; these are the data we report here. Two additional pointings
were centered 18′′ to the north and to the south of the first pointing; no
emission was detected at either offset pointing. In total we observed a field
of view of 16′′ $\times$ 56′′.

We reduced the data using the IDL version of the KCWI Data Extraction and
Reduction Pipeline (v1.2.1) and the IFSRED IDL library [1]. Sky subtraction
was performed using the IDL pipeline with manual selection of a sky mask
region within each pointing. The standard star BD+33d2642 was used for flux
calibration. Following the pipeline stages, we resampled the data onto 0.29′′
$\times$ 0.29′′ spaxel grids using the routine IFSR_KCWIRESAMPLE. The
resulting data cube has dimensions of 56 $\times$ 66 spaxels, covering 16.2′′
$\times$ 19.1′′. Assuming a $\Lambda$CDM Planck 2020 cosmology[2] the physical
dimensions of a single KCWI pointing correspond to 97 kpc $\times$ 114 kpc at
$z$=0.4512.

##  0.2 Spatially-integrated spectrum and line fluxes.

In addition to having spectra at each spaxel in the full datacube, we created
three spatially-integrated spectra by summing the observed flux in rectangular
apertures spanning the central 3$\times$3, 7$\times$7, and 15$\times$15
spaxels, corresponding to physical scales of 5, 12, and 26 kpc, respectively.
We detect strong [O ii] 3726Å and 3729Å and weak Mg ii 2796Å and 2803Å and [Ne
iii] 3869Å line emission in the spatially-integrated spectra (Extended Data
Fig. 1), which also include continuum emission from stars in the central
galaxy. The [O ii] doublet is blended in most spaxels, due to high velocity
dispersion. The [O ii] and [Ne iii] emission trace ionized gas; while Mg ii is
ionized gas, it is a resonant transition to the ground state and has an
ionization state that overlaps with neutral hydrogen, such that it traces
neutral gas[3]. Both [O ii] and Mg ii trace gas temperatures of $\sim 10^{4}$
K. The [O ii] doublet is observed at air wavelengths of 5407.2 and 5411.2Å;
the spectroscopic redshift of ORC4 is therefore $z=0.4512$. The [Ne v] 3426Å
emission line is not detected.

The spatially-integrated spectrum spanning the largest scales is used to
determine integrated fluxes. The flux of the [O ii] doublet, corrected for
Galactic extinction, is $8.8\ (\pm 0.1)\times 10^{-16}\ erg\ s^{-1}\ cm^{-2}$.
The total Mg ii flux is $0.64\ (\pm 0.13)\times 10^{-16}\ erg\ s^{-1}\
cm^{-2}$; the ratio of the [O ii] to Mg ii flux is 14. The restframe EW of the
Mg ii emission is 7 Å, and the observed ratio of the blue to red Mg ii
emission doublet is 1.9, close to the intrinsic flux ratio of 2.0, which
implies that the gas is optically thin[3] and not substantially impacted by
absorption. The flux of the [Ne iii] 3869 Å line is 0.72 ($\pm 0.11)\times
10^{-16}\ erg\ s^{-1}\ cm^{-2}$. The Mg ii and [Ne iii] emission are too weak
to detect in spatially resolved maps. In the spatially-integrated spectrum the
kinematics of the Mg ii and [Ne iii] emission are consistent with the [O ii]
emission, both in terms of velocity centroid and width.

##  0.3 Spectral cube analysis.

We fit Gaussian profiles to the observed emission lines, fitting Mg ii, [O
ii], and [Ne iii] in the spatially-integrated spectra and [O ii] in each
spaxel of the full datacube. As the [O ii] doublet is blended, we fix the [O
ii] flux ratio to be 1.2, corresponding to an electron density of $\sim 200\
cm^{-3}$. We did not tie the kinematics of the Mg ii, [O ii], and [Ne iii]
lines together. The continuum level around each line is determined locally
using the flux within 40 Å on either side of the emission singlet or doublet.

We define $V_{50}$ as the central velocity of the emission, as determined by
the centroid of the Gaussian fit. The velocity dispersion $\sigma$ is also
determined by the Gaussian fit, from which the instrumental resolution is
subtracted in quadrature. The ratio of $V_{50}$ to $\sigma$ in ORC4 is low,
with a median value of 0.19; this is an order of magnitude lower than what is
seen in star-forming galaxies[4]. The statistical errors on $\sigma$ range
from $\sim$10 km s-1 at the center of the nebula to $\sim$20 km s-1 at a
radius of 5 kpc and $\sim$65 km s-1 at a radius of 10 kpc. The systematic
errors on $\sigma$ are subdominant; changing the [O ii] flux ratio from
extreme values of 1.0 to 1.5 results in median differences to $\sigma$ of only
$\sim$2 km s-1 and does not change our results.

The systemic redshift of $z=0.4512$ is defined by the centroid of the observed
wavelengths of the [O ii] emission in the largest spatially-integrated
spectrum. The [O ii] EW is calculated by comparing the total flux in the [O
ii] doublet to the fit of the continuum level around [O ii], which is divided
by $(1+z)$ to obtain a restframe EW.

In Figs. 1 and 3 we show those spaxels for which the [O ii] flux has a signal
to noise ratio of 3 or greater. We calculate the surface brightness radial
profile using RADP in IDL and fit for the slope of the [O ii] profile from
$6-15$ kpc, beyond the region of the seeing disk. The errors shown in Fig. 2
are standard errors, calculated from the variation in the surface brightness
within that radial bin, and the errors in Fig. 3 are the standard errors
calculated from the variation in the central velocity fit across the central
five spaxels of the nebula.

##  0.4 Surface brightness profile and Sersic fit.

The central $\sim$5 kpc of the radial surface brightness profiles are not
resolved due to the observational seeing; from 6-15 kpc the best fit power law
slope to the [O ii] profile is $-2.8$. We fit the radial surface brightness
profile of both the [O ii] nebula and the stellar continuum with a Sersic
profile, convolving with the PSF of the observations. The Sersic profile is:

$I(r)=I_{e}\ exp\\{-(1.9992n-0.3271)\ [(r/r_{e})^{(1/n)}-1)]\\}$ (1)

where $I$ is the surface brightness as a function of radius, $I_{e}$ is the
surface brightess at the half light radius, $r_{e}$, and $n$ is the Sersic
index, which reflects the amount of curvature in the profile. We use CONVOL in
IDL to convolve this profile with a Gaussian with a FWHM of 0.9′′to account
for the seeing. The PSF-corrected best-fit Sersic index to the [O ii] emission
yields $n=1$, while for the stellar continuum the best fit is $n=2$. We
verified these results using GALFIT with a Gaussian PSF. The stellar emission
in the ORC4 central galaxy is more extended than the [O ii] emission; the
galaxy half light radius is 4.8 kpc.

##  0.5 Stellar mass and age estimation.

There is another galaxy near ORC4, as seen in Fig. 1. Long-slit spectroscopy
reveals that it is at the same redshift as ORC4 (D. S. N. Rupke et al., in
preparation) and is 36 kpc away; these galaxies are near enough to likely be
interacting[5]. We estimate stellar masses and ages for ORC4 and its companion
galaxy by fitting the spectral energy distributions (SEDs) using Prospector
[6], using a redshift of $z=0.4512$ for both galaxies (Extended Data Fig. 2).
Prospector utilizes gridless Bayesian parameter estimation and incorporates
Markov chain Monte Carlo (MCMC) posterior sampling with emcee to infer the
properties of the underlying stellar populations using Flexible Stellar
Population Synthesis models (FSPS)[7]. For ORC4, we fit broadband SDSS DR8[8]
ugriz and WISE[9] W1 and W2 photometry, and for the companion galaxy we fit
SDSS ugriz photometry. We assume a Kroupa initial mass function[10] and a
power law attenuation curve with index = -0.7.

For ORC4 we assume a parametric $\tau$-model star formation history with a
burst where the star formation rate is given as $SFR\sim\exp(t_{age}/\tau)$,
and we allow an AGN contribution. The free parameters in the model are the
galaxy stellar mass and age ($t_{age}$), the fraction of stellar mass formed
in the burst ($f_{burst}$), the fraction of the stellar age at which a
starburst occurred ($f_{age,burst}$), $\tau$, and the fraction of the light
that is due to an AGN. All timescales have units of Gyr. We set a log-normal
prior on the stellar mass with a mode of $\ln(M_{*}/M_{\odot})=23$ and
$\sigma=15$, as well as a log-uniform prior on $t_{age}$ such that
$0.1<t_{age}/\text{Gyr}<t_{universe}(z)$ where $t_{universe}(z)$ is the age of
the universe at the given redshift of $z=0.4512$. For the companion galaxy we
similarly assume a parametric $\tau$-model star formation history, but we do
not include a burst or an AGN contribution. Stellar mass and $t_{age}$ are the
only free parameters in the model for the companion galaxy, and we assume the
same priors that we implemented for ORC4. For both ORC4 and its companion
galaxy the stellar metallicity ($\log Z_{*}/Z_{\odot}$) is updated for each
random draw of $M_{*}$ within $\pm 0.1$ dex of the stellar mass-metallicity
relation[11].

The SED fit to the companion galaxy reveals that it is massive, with log
$(M_{*}/M_{\odot})=10.60\ (+0.10/-0.16)$, and old, with $t_{age}=7.0\
(+1.2/-1.4)$ Gyr, similar to what is found for the ORC4 central galaxy.

##  0.6 Radio luminosity conversion.

To convert the observed radio flux of the central source in ORC4 to a
luminosity at a restframe of 200 MHz, we use the following equation:

$L_{\nu}=4\pi
D_{L}^{2}(1+z)^{-(1+\alpha)}(\nu/\nu_{obs})^{(1+\alpha)}(\nu_{obs}/\nu)F_{\nu_{obs}}$
(2)

where $D_{L}$ is the luminosity distance, $\alpha$ is the observed spectral
index of ORC4 of $-0.92$, $\nu$ is the observed-frame frequency of the radio
luminosity we are deriving, here equal to 200 MHz $\times(1+z)$, $\nu_{obs}$
is the observed frequency of 325 MHz, and $F_{\nu_{o}bs}$ is the observed flux
of 1.43mJy. We use this same equation to convert to a restframe 1.4 GHz
luminosity.

##  0.7 Ionized gas mass estimate.

We use the [O ii] luminosity to estimate the mass of the warm ionized gas in
ORC4. Following Pandya et al. (2017)[2], we assume [O ii]/H$\beta=6$ to
compute an H$\beta$ luminosity and then use the following equation:

$M_{H\beta}=28.2\times 10^{8}L_{H\beta,43}\ n_{e,100}^{-1}\ M_{\odot}$ (3)

where $L_{H\beta,43}$ is the H$\beta$ luminosity in units of $10^{43}\ erg\
s^{-1}$ and $n_{e,100}$ is the electron density in units of 100 cm-3. We
assume $n_{e,100}$=1 to facilitate a comparison with the results of Pandya et
al. (2017). The ratio of [O ii]/H$\beta$ assumed is similar to the value found
in long-slit spectroscopy of ORC4 ($\sim$5.5, D. S. N. Rupke et al., in
preparation).

##  0.8 Alternative scenarios for ionized gas origin.

There are several alternative scenarios for the origin of the ionized gas in
ORC4. The first is that it is associated with the interstellar medium (ISM) of
the central galaxy. The velocity gradient in the ionized gas could be due to
rotation of the galaxy; many early-type galaxies have spatially extended,
rotating [O ii] emission[2]. However, the velocity gradient is asymmetric and
the radial extent and the velocity dispersion of the [O ii] emission in ORC4
are unusually high. From the [O ii] luminosity we estimate a warm ionized gas
mass of $3\times 10^{7}$ $M_{\odot}$, which is over an order of magnitude
higher than that found for massive early-type galaxies[2]. Additionally, the
observed [O ii] EW is far higher than what is seen in similar early-type
galaxies, such that it is very unlikely that this ionized gas in ORC4 is
simply cool gas in the ISM of the central galaxy.

Another possibility is that ORC4 may have a bipolar outflowing wind driven by
the central AGN, creating the extended [O ii] emission. “Red geysers” are
early-type galaxies that have bisymmetric EW maps with strong emission lines
that appear to trace large-scale outflows[12, 13]. Red geysers account for
5-10% of local, early-type galaxies[14]. AGN-driven outflows can result in
ionized gas velocity gradients, and the central galaxy in ORC4 does have an
AGN, such that it is plausible that it has an AGN-driven outflow. However,
such outflows are typically not circular and radio AGN outflows in particular
are jet driven and collimated[15]. Additionally, the [O ii] EW map of ORC4
does not display the collimated, bisymmetric observational signatures of red
geyser EW maps. Finally, the [O ii] luminosity of ORC4 is two magnitudes
higher than what is expected for a typical AGN with the same radio continuum
luminosity (Extended Data Fig. 3), and the SED fit to the optical and infrared
emission finds that only 3% of the light is from an AGN component. Thus
another mechanism is required to explain the strong, extended [O ii] emission,
and the high luminosity and EW in particular, in ORC4.

It is also unlikely that the [O ii] emission could result from a cooling flow
as ionized gas in cooling flows typically has a complex, and often
filamentary, morphology[16]. Additionally, in order for a cooling flow to be
sustained the circumgalactic medium should be well settled and not have any
disrupting events; the presence of the large-scale radio emission argues
against this. This galaxy also does not appear to be in a well-defined
overdensity or galaxy cluster [5] where cooling flows might be expected.

##  0.9 Outflowing galactic wind simulation.

We run the spherically-symmetric wind simulation using the Athena
hydrodynamics code. The galaxy is modeled as an isolated, spherical wind
injection region of radius $r=$ 7 kpc. The gravitational potential is provided
by a static NFW halo of virial mass 1013 $M_{\odot}$, which was chosen from
the low-redshift stellar mass to halo mass abundance matching relation [17]
for a galaxy of stellar mass $\sim 10^{11.2}$ $M_{\odot}$, similar to the
value estimated for ORC4’s central galaxy. The simulated circumgalactic medium
is initialized to be in isothermal hydrostatic equilibrium with $T=10^{6}$ K,
close to the virial temperature for a halo of mass $10^{13}$ $M_{\odot}$, and
with a gas density near the wind injection region of $10^{-26}\ g\ cm^{-3}$.
The circumgalactic medium is initialized with random velocity and density
perturbations with amplitudes on the order of 0.2 times the virial velocity
and 0.3 times the density, respectively. The wind turns on immediately at the
beginning of the simulation and blows spherically-symmetrically with a
velocity of 450 km s-1 and with a mass outflow rate of 200 $M_{\odot}$yr-1,
and then the wind turns off at t = 200 Myr. We run just one octant of the
sphere to reduce computational costs. The simulation domain is a cube with
side length 350 kpc where the wind injection region is located in one corner,
at the origin. There are three levels of static mesh refinement: within a 44
kpc cube cornered at the origin, and thus closest to the wind injection
region, the spatial resolution is 438 pc. The resolution decreases to 877 pc
out to the edges of a 175 kpc cube, and then again to 1.75 kpc out to the edge
of the simulation domain at 350 kpc.

In the simulation that most closely matches the observational inputs, at 740
Myr after the starburst begins, or 540 Myr since the star formation and
resulting wind have ceased (as the wind blows for 200 Myr), the forward shock
is located at 215 kpc in radius and has a Mach number of 1.3, consistent with
the shock Mach number required to produce the observed synchrotron emission in
ORCs[19]. At this same time, the previously-shocked wind that gives rise to
the [O ii] emission is falling back toward the galaxy and extends from a
radius of 0 to 35 kpc. The velocity dispersion of the previously-shocked wind
gas that is returning to the galaxy is $\sim$120 km s-1 within a 20 kpc
radius, somewhat lower than the [O ii] velocity dispersion observed in ORC4.
The velocity dispersion in this simulation is a lower limit on the expected
velocity dispersion of [O ii] emitting gas, as the simulation is not resolving
the chaotic series of shocks that are likely producing the [O ii] emission
within the simulated $10^{4}$ K shocked wind region. The parameters of the
wind (velocity, mass outflow rate, the timescale that the wind is “on”) and
the halo (virial mass, temperature, density) were chosen to produce results
that match the properties of ORC4, though we found that a wide range of
parameters were generically able to produce a forward shock on large scales
and shocked wind falling back to the galaxy on small scales that could be
observed hundreds of Myr after the wind from the galaxy had ceased. The
simulation does not predict observed [O ii] emission line luminosities or the
ionized gas mass in the shocked wind gas due to the computational difficulty
of radiative transfer modeling.

##  0.10 Data and Code Availability.

Raw data generated at the Keck Observatory are available at the Keck
Observatory Archive (koa.ipac.caltech.edu) following the standard 18-month
proprietary period after the date of observation. The reduced KCWI spectral
datacube and the results of the [O ii] emission line fits are available on
Zenodo at https://doi.org/10.5281/zenodo.8377942. The code used is available
upon request from the first author.

## References

  * [1] Rupke, D. S. N.  IFSRED: Data Reduction for Integral Field Spectrographs.  Astrophysics Source Code Library, ascl:1409.004 (2014).
  * [2] Planck Collaboration, Aghanim, N. _et al._ Planck 2018 results. VI. Cosmological parameters.  _A &A_ 641, A6 (2020).
  * [3] Chisholm, J., Prochaska, J. X., Schaerer, D., Gazagnes, S. & Henry, A.  Optically thin spatially resolved Mg II emission maps the escape of ionizing photons.  _MNRAS_ 498, 2554–2574 (2020).
  * [4] Green, A. W. _et al._ DYNAMO - I. A sample of H$\alpha$-luminous galaxies with resolved kinematics.  _MNRAS_ 437, 1070–1095 (2014).
  * [5] Norris, R. P., Crawford, E. & Macgregor, P.  Odd Radio Circles and Their Environment.  _Galaxies_ 9, 83 (2021).
  * [6] Johnson, B. D., Leja, J., Conroy, C. & Speagle, J. S.  Stellar Population Inference with Prospector.  _ApJS_ 254, 22 (2021).
  * [7] Conroy, C., Gunn, J. E. & White, M.  The Propagation of Uncertainties in Stellar Population Synthesis Modeling. I. The Relevance of Uncertain Aspects of Stellar Evolution and the Initial Mass Function to the Derived Physical Properties of Galaxies.  _ApJ_ 699, 486 (2009).
  * [8] Aihara, H. _et al._ The Eighth Data Release of the Sloan Digital Sky Survey: First Data from SDSS-III.  _ApJS_ 193, 29 (2011).
  * [9] Wright, E. L. _et al._ The Wide-field Infrared Survey Explorer (WISE): Mission Description and Initial On-orbit Performance.  _AJ_ 140, 1868–1881 (2010).
  * [10] Kroupa, P.  On the variation of the initial mass function.  _MNRAS_ 322, 231 (2001).
  * [11] Ma, X. _et al._ The origin and evolution of the galaxy mass-metallicity relation.  _MNRAS_ 456, 2140 (2016).
  * [12] Cheung, E. _et al._ Suppressing star formation in quiescent galaxies with supermassive black hole winds.  _Nature_ 533, 504 (2016).
  * [13] Roy, N. _et al._ Evidence of Wind Signatures in the Gas Velocity Profiles of Red Geysers.  _ApJ_ 913, 33 (2021).
  * [14] Roy, N. _et al._ Detecting Radio AGN Signatures in Red Geysers.  _ApJ_ 869, 117 (2018).
  * [15] Blandford, R., Meier, D. & Readhead, A.  Relativistic Jets from Active Galactic Nuclei.  _ARA &A_ 57, 467–509 (2019).
  * [16] McDonald, M., Veilleux, S., Rupke, D. S. N. & Mushotzky, R.  On the Origin of the Extended H$\alpha$ Filaments in Cooling Flow Clusters.  _ApJ_ 721, 1262–1283 (2010).
  * [17] Girelli, G. _et al._ The stellar-to-halo mass relation over the past 12 Gyr. I. Standard $\Lambda$CDM model.  _A &A_ 634, A135 (2020).
  * [18] Best, P. N. & Heckman, T. M.  On the fundamental dichotomy in the local radio-AGN population: accretion, evolution and host galaxy properties.  _MNRAS_ 421, 1569 (2012).

We thank Fabrizio Brighenti for useful feedback on an earlier draft of the
paper, and we thank Ray Norris and Huib Intema for sharing the GMRT discovery
image of ORC4. A.L.C. acknowledges support from the Ingrid and Joseph W.
Hibben chair at UC San Diego. C.L. thanks Charles Cimino IV for computing
resources and assistance. The data presented herein were obtained at the W. M.
Keck Observatory, which is operated as a scientific partnership among the
California Institute of Technology, the University of California, and the
National Aeronautics and Space Administration. The Observatory was made
possible by the financial support of the W. M. Keck Foundation. The authors
wish to recognize and acknowledge the very significant cultural role and
reverence that the summit of Maunakea has always had within the indigenous
Hawaiian community. We are most fortunate to have the opportunity to conduct
observations from this mountain.

D.S.N.R. and A.L.C. conceived the observations, following C.A.T.’s suggestion
that ORCs may be similar to the galaxies we have been studying. A.L.C.
obtained the observing time, S.P. performed the KCWI observations and led the
data reduction. A.L.C. performed all of the data analysis, with input from
D.S.N.R. C.L. ran the starburst-wind model simulation, and K.W. performed the
SED fitting. A.L.C. wrote the manuscript, with input from all coauthors.
Figures were created by S.P., A.L.C., K.W., C.L., and R.H.

Figure 1: Spatially-integrated spectrum of the gas nebula and host galaxy
stellar continuum emission in ORC4. (a) Strong [O ii] 3726, 3729 Å and weak Mg
ii 2796, 2803 Å and [Ne iii] 3869 Å emission are observed in a spatially-
integrated spectrum of ORC4 (spanning the inner 26x26 kpc), at wavelengths
that correspond to $z=0.4512$; weak stellar continuum is also detected. The
solid black line shows the observed spectrum, the pink dotted line is the
1$\sigma$ error spectrum, and the orange dashed vertical lines indicate the
observed wavelengths of Mg ii, [O ii], and [Ne iii] at the spectroscopic
redshift of the source. (b) Gaussian fits to the Mg ii and [O ii] emission
doublets and the [Ne iii] singlet in the spatially-integrated spectrum. Blue
and cyan lines show fits to the individual emission lines in each doublet.
Figure 2: Spectral energy distribution (SED) fit to optical SDSS and near
infrared WISE photometry of ORC4. (a) The observed photometry of ORC4 is shown
as green circles with 1 $\sigma$ error bars, while orange circles show the
photometry implied from the best fit stellar population model, including an
AGN contribution, which is shown in blue. Flux values are given in mJy and
observed frame wavelengths in $\mu$m. (b) The derived distributions from the
SED fit for the stellar mass, stellar age, star formation history parameters,
and AGN contribution for ORC4, as well as the covariance between the
parameters. Figure 3: Comparison of ORC4 radio continuum and [O ii] line
luminosity to radio AGN. The 1.4 GHz radio continuum luminosity and [O ii]
emission line luminosity of ORC4, shown with a red circle, compared to the
radio-loud AGN sample of Best & Heckman (2012)[18], shown with grey points.
The [O ii] luminosity of ORC4 is two orders of magnitude higher than the
median value for an AGN with the same 1.4 GHz radio continuum luminosity.


# Single atom chemical identification of TMD defects in ambient conditions

Dunn E. J.1,2, Young R. J.1 and Jarvis S. P.1

(1Physics Department, Lancaster University, Lancaster, LA1 4YB, UK.
2Univ. Lille, CNRS, Centrale Lille, Univ. Polytechnique Hauts-de-France, UMR
8520 - IEMN - Institut d’Electronique de Microélectronique et de
Nanotechnologie, F-59000 Lille, France.
)

###### Abstract

The presence of defects in transition metal dichalcogenides (TMDs) can lead to
dramatic local changes in their properties which are of interest for a range
of technologies including quantum security devices, hydrogen production, and
energy storage. It is therefore essential to be able to study these materials
in their native environments, including ambient conditions. Here we report
single atom resolution imaging of atomic defects in MoS2, WSe2 and WS2
monolayers carried out in ambient conditions using conductive atomic force
microscopy (C-AFM). By comparing measurements from a range of TMDs we use
C-AFM to chemically identify the most likely atomic species for the defects
observed and quantify their prevalence on each material, identifying oxygen
chalcogen substitutions and transition metal substitutions as the most likely,
and most common, defect types. Moreover, we demonstrate that C-AFM operated in
ambient environments can resolve subtle changes in electronic structure with
atomic resolution, which we apply to WSe2 monolayers doped using a nitrogen
plasma, demonstrating the capability of C-AFM to resolve electronic, and
chemical-specific, details at the atomic scale.

## Introduction

The properties of transition metal dichalcogenides (TMDs) make them promising
candidates for improving existing technologies and enabling those of the
future [1, 2]. TMDs have a range of interesting properties and applications;
some monolayers are direct band-gap semiconductors with excellent optical
properties [3], they can be highly catalytic as a result of chalcogen defects
introducing active sites[4, 5, 6] and, as layered materials, they can be
intercalated, allowing the storage of ions in batteries [7], or stacked to
form heterostructures with new properties [8]. The disorder introduced by
defects can also be harnessed for the creation of security devices [9],
exploiting the presence of the unique arrangement of defects in individual TMD
flakes which is reflected in the photoluminescence response [10, 11].

Scanning probe techniques offer a means of direct real-space imaging providing
insight into the defects present after fabrication and the processes through
which they are created and transformed. Recently the ability to resolve these
defects with atomic resolution has been demonstrated on TMDs and MXenes using
conductive atomic force microscopy (C-AFM) [12, 13] and higher eigenmode
techniques [14], complementing the traditional Scanning Tunnelling Microscopy
(STM) [15, 16, 17, 18, 19, 20, 21, 22, 23].

In this work, the intrinsic defects of mechanically exfoliated TMD monolayers
are characterised using C-AFM with atomic resolution in ambient conditions. By
comparing results from MoS2, WSe2 and WS2 monolayers, we are able to identify
the most likely chemical origins for single atom defects and quantify their
prevalence on each material. Moreover, we demonstrate that C-AFM operated in
ambient environments can not only resolve defects with atomic resolution, but
also detect subtle changes in electronic structure. We apply this method to
WSe2 monolayers intentionally doped through exposure to nitrogen plasma, where
we show that substitutions of individual oxygen and nitrogen atoms on
chalcogen sites (OSe and NSe) can be distinguished.

## Methods

## Sample preparation.

To fabricate TMD monolayer samples bulk TMD crystals (MoS2 \- Manchester
Nanomaterials, WSe2 \- HQ Graphene, and WS2 \- 2D Semiconductors) were
mechanically exfoliated onto Nitto Blue Tapes (Acrylic/PVC) where repeated
exfoliation dispersed and thinned the flakes. The flakes were then exfoliated
from the tapes onto pre-prepared PDMS stamps affixed to glass slides where
flakes were then assessed for quality using an optical microscope. Few-layer
flakes were selected for further study with Raman and photoluminescence (PL)
in a Horiba labRAM instrument using a Laser Quantum 532 nm laser. Spectra were
used to identify monolayers which were then transferred using the PDMS dry
transfer technique onto a freshly exfoliated highly ordered pyrolytic graphite
(HOPG) substrate.

Here we use HOPG which exhibits a near ohmic density of states and can be
prepared with large atomically flat terraces [24]. In particular, good
electrical contact with HOPG ensures that gap states within the TMD materials
can be easily measured with C-AFM.

Plasma treated samples were prepared in a custom-built glass barrelled plasma
reactor. The glass barrel of the reactor was a QVF process pipe (De Dietrich,
UK) – dimensions $500$ mm length, $100$ mm diameter, clamped between two
custom-made steel end plates which serve as ground electrodes. The system was
evacuated using an Edwards RV3 rotary vane pump (base pressure: c. $2x10^{-3}$
mbar). The pressure within the system was monitored by a Thermovac TTR 96N SC
Pirani gauge and Display One Controller (Leybold UK Ltd) and an in-line liquid
nitrogen cold trap was installed to enhance the system base pressure. Plasma
was ignited using a RFG-C-100-13 power generator and automatic matching
network (Coaxial Power System Ltd, Eastbourne, UK), operated at a frequency of
$13.56$ MHz, connected to a $1$ cm thick copper braid (Tranect Ltd, Liverpool,
UK) wrapped around the glass barrel three times.

## Conductive atomic force microscopy.

All scanning probe microscopy (SPM) measurements were collected with a Bruker
MultiMode 8 equipped with a Nanoscope V controller housed in the Lancaster
ultra-low-noise IsoLab facility. The IsoLab facility is comprised of a
dedicated building with three above-ground laboratories each contained in
their own separate pod. In the basement of each pod sits a 50-ton concrete
isolation block with passive isolation providing exceptional low-vibration
environments. For the flattening of monolayer flakes and removal of
contaminants, the AFM was operated in contact mode using NuNano Scout 70
probes ($f_{0}$ $\sim 70$ kHz; nominal spring constant 2 Nm-1). Samples were
typically then checked in PeakForce mode to ensure no damage to flakes.
Conductive atomic force microscopy (C-AFM) was carried out using Nanoworld
CDT-FMR conductive diamond-coated Si probes ($f_{0}$ $\sim 105$ kHz; nominal
spring constant 6.2 Nm-1). C-AFM measurements were collected by operating in
contact mode at a deflection setpoint of 5 nm, corresponding to a force of
$\sim 31$ nN. A bias voltage was applied to the sample, with the resulting
current measured via the conductive probe using a FEMTO DLPCA-200 amplifier.

## Results.

## Flattening of monolayer flakes.

Atomic resolution imaging with conductive atomic force microscopy (C-AFM)
requires atomically flat and electrically conductive substrates, without which
features cannot be resolved.

Following transfer, the TMD monolayer flakes were flattened using the so-
called ‘Nano-Squeegee’ technique [25]. This process removes the majority of
contamination underneath the TMD flakes introduced from the transfer method,
greatly improving the electrical contact and creating atomically flat regions
of TMD for high-resolution C-AFM measurements. Flattening is achieved by
operating the AFM in contact mode at a higher deflection set point and
repeatedly scanning the surface, ensuring that the TMD flakes remain undamaged
(see Figure S1).

## Atomic defects in MoS2.

The formation and structure of atomic defects on MoS2 is very well studied
[12, 17, 20, 22, 23, 26, 27, 28, 29], and can be used to aid assignment of the
defects observed with C-AFM in other materials. In Figure 1 we show the three
main defects observed with C-AFM carried out in ambient conditions including
images acquired with true atomic resolution.

Figure 1:  Single atom defects in MoS2 measured in ambient conditions. (a)
C-AFM image showing atomic scale defects visible as either bright or dark
features located on sulfur lattice sites. (b) High resolution images show
clear augmentation of the measured current in the nearest neighbour sites
around individual defects, with line profile analysis indicating an
augmentation of almost $50$%. (c) Ball-and-stick cartoons depicting the most
common defects in MoS2. (d) C-AFM image showing ‘diffuse’ defects with large
radius. All C-AFM images in this figure were collected on a MoS2 monolayer on
HOPG collected at +$0.5\,$V sample bias.

In Figure 1a) we show the two most common defects observed, identified as
either bright or dark features one atom in size. Furthermore, the surrounding
nearest neighbour atoms for these defects appear brighter (increased current)
compared to the rest of the MoS2 flake. This is most clearly observed for the
dark point defects, as shown in more detail in Figure 1b). The line profile in
Figure 1b) reveals that the current signal is almost 50% larger for the
nearest neighbour atoms compared to atoms further from the defect, suggesting
an increase in local unoccupied states due to the presence of the atomic
defect, consistent with previous observations on, and simulations of, TMDs
[15, 17, 20, 23, 30]. The precise origin of these single-atom size defects is
not universally agreed, with reports suggesting that either sulfur vacancies
[17, 12, 13, 26] or oxygen substitutions [15, 22, 31, 32] are responsible.
Zhang et al. [33] suggest that their most common defects are caused by
transition metal vacancies because they observe the defect to lie in-between
chalcogen sites, but this is inconsistent with our data shown in Figure 1. The
data from Barja et al. is particularly compelling, with the combination of
STM, dI/dV, and nc-AFM strongly supporting oxygen substitution. Furthermore,
the results from Barja et al. are similar to our own assignment that the
bright and dark atomic defects arise from sulfur defects in the bottom and top
layers of the MoS2 monolayer, respectively, with enhancement in the bright
defect arising from stronger interaction with the substrate. We also note that
whilst our samples are stored under vacuum, ambient measurements necessitate
significant exposure to atmosphere and so atomic vacancies are unlikely.
Several studies suggest that oxygen substitutions are energetically favourable
[34, 35] and that the presence of a chalcogen vacancy encourages the
adsorption and dissociation of molecular oxygen; so that in the presence of
oxygen, chalcogen vacancies assist in their own transformation into oxygen
substitutions [27, 28, 35].

Another less commonly observed defect is shown in Figure 1d). These more
diffuse defects appear as either bright (Figure 1d) or dark (Figure S2)
features and are characterised by larger lateral dimensions extending over
multiple lattice sites. Diffuse defects have previously been observed in C-AFM
and attributed to transition-metal substitutions [12, 26]. Similar dark $5$ nm
radius defects were imaged via STM by Addou et al. who also see significant
variation between samples that have been prepared in the same way [17].
Diffuse defects have also been observed in STM by Schuler et al. and
attributed to charged defects [36]. Charged defects induce band bending,
however, in this case the shift is always in the same direction for both the
conduction and valence bands. This means that an augmentation of the current
at positive (negative) bias leads to a depletion of the current at negative
(positive) bias. Figure 2b-e) shows that the diffuse defects we observe
consistently deplete the current from $-1$V to $1$V, therefore we rule out
charged defects as the origin of this effect.

Transition metals can have similar or favourable binding energies when
substituted at a molybdenum site when compared to the pristine lattice [37].
Ta, W, Nb, Hf and Zr have higher binding energies than Mo and, as a result, if
they are present they could be expected to replace the Mo. We find that
diffuse defects are not observed as commonly on MoS2 as the single atom size
defects. This is consistent with the expectation that transition-metal
substitutions are relatively rare, occurring only during initial growth of the
bulk MoS2 crystal. These observations therefore support an assignment of
diffuse defects to transition metal substitutions, similar to previous reports
[12, 26].

## Atomic defects in WSe2.

Beyond MoS2, two of the most commonly studied 2D-TMDs are WSe2 and WS2. There
are many STM studies into the chemical nature of defects on WSe2 and WS2 in
synthetic samples [13, 15, 18, 19, 21, 31, 36, 38, 39]. Here we apply C-AFM as
described above in order to identify the frequency and type of atomic vacancy
observed in mechanically exfoliated monolayers of each material. Figure 2a)
shows a typical atomic resolution C-AFM image collected on WSe2 where single
atom vacancies are observed. In this case, we observe a much greater number of
atomic defects as compared to MoS2, most of which appear as dark features
which are attributed to Se vacancies or OSe substitutions as described in the
previous section. The observation of predominantly dark single atom vacancies
is consistent across measurements, implying a greater number of defects in the
top layer of the WSe2. We tentatively attribute this to increased atmospheric
exposure in the top layer, resulting in a larger number of OSe substitutions.

Figure 2:  Atomic resolution C-AFM measurements of defects in WSe2 and WS2.
(a) C-AFM image of atomic depleting defects on a WSe2 monolayer imaged at
$1\,$V. C-AFM image of a diffuse defect on a WSe2 monolayer imaged at (b)
$-1\,$V and (c) $1\,$V sample bias. (d) A profile across the diffuse defect
indicates a diameter of approximately $5$ nm. C-AFM images on WS2 are
dominated by a Moiré pattern from the lattice mismatch with the underlying
HOPG. (e) Augmenting (white circle) and depleting (dashed white circle)
defects on a WS2 monolayer observed at $0.5$ V sample bias. Images collected
at (f) $0.5$ V and (g) $-1$ V sample bias in another region show only the
Moiré pattern, with no visible defects. (h) An example of defect, or
potentially molecular adsorbate, visible at $-0.8$ V sample bias.

The increased number of single atom features observed on WSe2 appears to
support reports that defects are more numerous in WSe2 due to the reduction in
air stability for chalcogens with higher atomic numbers [27, 40, 41]. The
comparison between MoS2 and WSe2 also allows us to provide more insight into
the assignment of the atomic defects as oxygen substitutions or chalcogen
vacancies. Simulated studies of MoS2 and WSe2 reveal similar formation
energies for chalcogen vacancies in each material. This is due to an increase
in binding energy when changing from Mo to W, which acts to counteract the
decrease in binding energy when S is exchanged for Se [27]. Comparatively, the
calculated energy barrier for oxygen adsorption onto WSe2 is much less than
that required for MoS2 [28], suggesting a preference for oxidation on WSe2.
The increased defect density we observe on WSe2 compared to MoS2, therefore,
suggests that oxidation is most likely responsible, further supporting the
assignment of oxygen substitutions rather than Se vacancies.

A higher density of diffuse defects is also observed on WSe2. As with other
samples, diffuse defects are still much less common than site-specific
defects; an observation in agreement with Yankowitz at al. [19]. A pair of
C-AFM images collected at $\pm 1$ V bias show one such defect in Figure 2b-c).
These defects have a similar lateral size to those on MoS2 and, similarly,
they may be found as both dark or bright features (Figure S2); though an
individual defect is consistently dark or bright across the voltage range
used. We note that all C-AFM images presented in this paper use a colour scale
that depicts larger current magnitudes as bright, and smaller current
magnitudes as dark. It is unclear why some diffuse defects deplete the density
of states rather than augment it, though many different transition metals
could potentially be substituted into the lattice.

In Figure 2c) which was collected with a sample bias of $1\,$V we also observe
defects with smaller radii. These defects are similar to the single atom size
defects observed in other images but appear to be around 1-2 nm in size. We
believe this is the result of a slightly blunt tip used in Figure 2b) and c),
as later on in the experiment the resolution improved, allowing us to see
these defects.

## Atomic defects in WS2.

Figure 2(e-h) shows typical C-AFM images collected on WS2 monolayers also
prepared on HOPG. In this case, strong Moiré patterns between the WS2 and the
HOPG substrate dominate C-AFM images. Repeat measurements suggest that Moiré
patterns for WS2:HOPG are almost ubiquitous, and are very rarely absent. This
can be compared to WSe2:HOPG, where Moiré patterns were also observed, but
only at particular bias voltages that could be easily avoided (Figure S3).

Despite complications arising from the Moiré patterns, defects are still
clearly observed, albeit in much smaller numbers compared to the other TMD
materials studied. Figure 2(e-h) shows high resolution C-AFM images collected
at a range of sample bias voltages between -1 V and 0.5 V. Small site specific
defects are highlighted in Figure 2(e) showing an augmenting (cyan circle) and
depleting defect (red circle) in a C-AFM image collected at 0.5 V sample bias.
In this case the Moiré pattern prevents any possibility for atomic resolution,
and the defects instead appear to influence a lateral region of $\sim 1$ nm.
Figure 2(g,h) shows two images from a separate measurement session with only
the Moiré pattern visible collected at 0.5 V and -1.0 V sample bias
respectively. A further example of the difficulty in identifying defects can
be found in Figure 2(h) which shows an example of another site specific defect
at -0.8 V sample bias. This defect appears to change lattice site during
scanning (slow scan direction is top-to-bottom), showing a ‘slice’ effect
where the defect switches from a dark to a bright feature, laterally offset by
around 0.5 nm. This leads to the conclusion that the defect in Figure 2(f) may
be a surface adsorbate, and so is shown here only for completeness. Diffuse
defects, as discussed for WSe2 and MoS2, were not observed on WS2. We
tentatively ascribe this to the difficulty in identifying such defects against
the background of the Moiré pattern, and, as discussed below, the relatively
low density of defects observed on WS2 compared to the other TMDs studied.

## Quantification of defect densities.

In Figure 3 WSe2, MoS2 and WS2 are compared side by side. Over a $30$ nm by
$30$ nm area the defect densities are $1.4\pm 0.1\times 10^{13}$, $1.6\pm
0.4\times 10^{12}$ and $6\pm 3\times 10^{11}$ defects per cm2 for WSe2, MoS2
and WS2 respectively. The densities of defects observed in Figure 3 are
consistent with expectations from density functional theory and molecular
dynamics studies of oxygen adsorption and dissociation on TMD surfaces and the
formation energy of vacancies. Compared to MoS2, WS2 should exhibit fewer
defects due to an increase in the vacancy formation energy arising from the
change of Mo to W [27]; which then, in turn, decreases the adsorption of
oxygen onto the surface [35]. Similarly, changing chalcogen from S in WS2 to
Se in WSe2 increases the vacancy formation energy with the same results.
However, it also leads to a higher energy barrier for oxygen adsorption [28].

Figure 3:  Quantification of defects on MoS2, WSe2 and WS2. $30$ nm by $30$ nm
current channel C-AFM images of site specific defects on (a) MoS2 collected at
$0.5$ V sample bias, (c) WSe2 collected at $1$ V sample bias, and (d) WS2
collected at $0.5$ V sample bias. (b) Shows a section of (a) in greater
detail. The image on WSe2 displays the most defects (125), followed by MoS2
(16) and then WS2 (6), however the Moiré pattern may hide some of the defects
in WS2. In all three materials a mix of augmenting and depleting defects are
observed

## Nitrogen doping of WSe2 via plasma.

The atomic resolution capability of C-AFM, and moreover, the sensitivity to
minor changes in conduction, suggests that C-AFM should be capable of
identifying variations in local conductance by dopant chemistry. WSe2
monolayers on HOPG were therefore introduced into a N2 plasma, creating NSe
substitutions that can be compared with the pre-existing OSe discussed above.
XPS measurements have previously confirmed that Mo-N bonds can be introduced
to MoS2 through plasma exposure, implying that N2 plasma creates nitrogen
substitutions on the chalcogen sites [29, 42], which have a similar formation
energy to OSe [43].

Samples were prepared using a bespoke plasma reactor designed for plasma
polymerisation experiments which require cleaner vacuum conditions than a
typical plasma oven. The number of NSe substitutions can be carefully
controlled by varying gas pressure, power, and the duration of plasma
treatment. Using conditions of 10 Pa pressure and 10 W power at 13.56 MHz, we
found that 10 s of exposure produced a similar density of NSe defects compared
to the pre-existing OSe. Figure 4a shows a C-AFM image of single atom OSe
depleting defects on WSe2:HOPG observed before plasma exposure. Following 10 s
exposure, a second type of defect appears as shown in Figure 4b, which once
again appears as a single atom depleting site, however, compared to OSe, a
faint atomic feature is still visible. Profiles taken across the defects are
shown in Figure 4(e) and (f) alongside enlarged C-AFM images of OSe and NSe in
Figure 4(c) and (d), respectively. The line profiles clearly show an increased
signal in the depleting region for NSe compared with OSe, which we find is
consistent across samples (and increase in number with increased plasma
exposure). In addition, we note that prior to plasma treatment, defects were
not frequently observed at nearest-neighbour chalcogen sites, however, in
Figure 4b paired defects are more common and can be observed for both ‘faint’
and single atom depleting defects. Enlarged examples are shown in Figure 4g-h.

It has been reported that nitrogen substitutions should add in-gap states near
the valence band [43], however, this doesn’t appear to be reflected in our
C-AFM measurements. We postulate this is because unlike OSe, NSe is not iso-
electronic, instead acting like a p-type dopant. The underlying HOPG provides
access to a carrier reservoir that can fill states, resulting in additional
electrons at the dopant site and a local negative charge. As previously
discussed, a key signature of a charged defect is asymmetry between negative
and positive bias, therefore, at positive bias a negatively charged defect
should be observed as a dark feature, which is consistent with the reduction
in current characterising the ‘faint’ defect sites at 1 V. This idea is
supported by similarities to charged defects observed in STM studies of CVD
TMDs with CHSe and CHS substitutions which are thought to be similar to nature
to NSe substitutions[44].

Figure 4:  Identification of defects in N2 plasma treated WSe2. Current
channel C-AFM images of WSe2 monolayer on HOPG imaged at $1$ V (a) before and
(b) after plasma exposure. Enlarged C-AFM images show (c) OSe and (d)
NSefeatures in detail, alongside line profiles shown in (e) and (f),
respectively. (g) C-AFM image of a pair of OSe defects, compared to (h) a
C-AFM image of a pair of NSe defects. Arrows are used throughout the figure to
highlight OSe (red) and NSe (blue) features.

## Mechanism of atomic resolution C-AFM in ambient conditions.

The diamond-coated probes used in this study last long enough for a full day
of measurements, whereas Pt-coated probes could only acquire a small number of
images before the conductive signal was lost. Therefore for practical reasons
diamond-coated probes were used. Initially, atomic resolution was not
expected, however, the diamond probes used achieved ‘true’ atomic resolution
in approximately $40$% of C-AFM imaging sessions, and lattice resolution
reliably achieved over $90$% of the time. We note that the probes remained
sufficiently sharp for atomic resolution despite a colossal loading force
compared to other high resolution imaging techniques [45]. We expect that the
high stability of the polycrystalline diamond coating plays an important role
in maintaining the sharpness of the tip apex. As discussed above, that atomic
resolution is not always achieved and despite the high stability of the
diamond probes ‘tip changes’ are often observed; generally resulting in a
change in image quality and potentially atomic resolution. Tip changes are
observed as sudden, discontinuous jumps in the current which originate from
changes in the nanostructures at the tip apex or from picking up and dropping
molecules on the surface. In such cases where atomic resolution is not
achieved, we find that the majority of defects appear as bright features
around 1-2 nm in size, and we no longer observe a mixture of bright and dark
point defects. This appears to be consistent with previous reports on CVD WS2
using C-AFM with nanoscale resolution [38]. We also note that whilst attempts
were made to measure I(V) at specific defect sites, this data was not reliable
due to the high thermal drift in ambient conditions.

High-resolution AFM in ambient conditions has been reported by a number of
groups[12, 13, 14, 26, 46, 47, 48, 49, 50, 51] with various proposed imaging
mechanisms. On TMD samples both Sumaiya et al. [12] and Nowakowski et al. [13]
attribute their atomic resolution images to an insulating layer around the tip
which is broken in a sufficiently small area as to create a single atom
pathway into the tip. Nowakowski et al. cite studies of the size of the
contact area relative to the area of the conducting channel by Enachescu et
al. [52] and Celano et al. [53]. Enachescu et al. tentatively propose an
insulating layer and Celano et al. show that the conductive area of the tip is
smaller than the contact area when tunnelling through SiO2 at tip loads from
$\sim 50$-$150$ nN. We note that all of the tip loads used by the reports
where true atomic resolution has been realised all lie below this range: $\sim
2$-$2.5$ nN[13], $\sim 0$ nN [12] and $\sim 31$ nN (this work). We also note
that atomic resolution in UHV-AFM experiments using Si probes have been
extensively studied [54, 45, 55, 56, 57], where it has long been known that a
single atom at the tip apex determines the tip-sample interaction [58] (and
therefore, what one might consider the tip-sample ‘contact’ in C-AFM), with no
need for insulating layers. We therefore propose that an insulating layer is
not necessary to achieve atomic resolution in ambient, and that the diamond
probe likely contains nanocrystallites at the tip-apex, which behave in a
similar way to the nano-asperities known to facilitate atomic resolution in
UHV-AFM experiments. Indeed, the type of diamond-coated probes used in this
study have been observed to terminate in nanoscrystallites in the $10$ nm
regime (See Figure S4).

## Conclusions.

We have shown that C-AFM is capable of routinely collecting atomic resolution
images of a variety of TMD monolayers prepared and measured in ambient
conditions. Not only is C-AFM sensitive to single atom defects, but also
subtle changes in the local electronic structure including that of nearest
neighbour atoms. MoS2 mechanically exfoiliated monolayers on HOPG substrates
are used as a well-studied reference where we identify OS and transition metal
substituents as the most likely, and most common defect types. Comparison
across MoS2, WSe2, and WS2, monolayers reveal consistent variation in defect
density of $1.4\times 10^{13}$, $1.6\times 10^{12}$ and $6.7\times 10^{11}$
defects per cm2 for WSe2, MoS2 and WS2 respectively, which correspond well to
expected formation energies. Finally, we show that minor variations in
electronic structure due to OSe and NSe substitutions can be measured in N2
plasma treated WSe2:HOPG monolayers further demonstrating the capability of
C-AFM to resolve electronic, and indeed chemical, details at the atomic scale.

By introducing elements into TMDs intentionally via methods such as the plasma
exposure demonstrated here it will be possible to find further evidence to
illuminate the chemical origin of the defects observed in this work, identify
non-intrinsic defects and establish the impact that these defects have on
physical properties of TMDs including the optical emission and catalytic
behaviour.

## Acknowledgements.

E.J.D. acknowledges support and funding from the GrapheneNOWNANO Centre for
Doctoral Training, EP/L01548X/1. R.J.Y. thanks The Royal Society grant number
UF160721, Air Force Office of Scientific Research grant number
FA9550-19-1-0397 and Engineering and Physical Sciences Research Council, Grant
No. EP/K50421X/1. S.P.J. thanks The Royal Society for grant PI70026, and the
Engineering and Physical Sciences Research Council (EPSRC) for grants
EP/V00767X/1 and EP/X026876/1, Quantum engineering of energy-efficient
molecular materials, QMol.

## References

  * [1] Manzeli S, Ovchinnikov D, Pasquier D, Yazyev OV, Kis A.  2D Transmission Metal Dichalcogenides.  Nature Reviews Materials. 2017;2:1–15.
  * [2] Autere A, Jussila H, Dai Y, Wang Y, Lipsanen H, Sun Z.  Nonlinear Optics with 2D Layered Materials.  Advanced Materials. 2018;30(1705963).
  * [3] Mak KF, Lee C, Hone J, Shan J, Heinz TF.  Atomically Thin MoS 2: A New Direct-Gap Semiconductor.  Physical Review Letters. 2010;105(136805).
  * [4] Zhu J, Yang R, Zhang G.  Atomically thin transition metal dichalcogenides for the hydrogen evolution reaction.  ChemMaterPhys. 2022;1:102–111.
  * [5] Ye G, Gong Y, Lin J, Li B, He Y, Pantelides ST, et al.  Defects Engineered Monolayer MoS2 for Improved Hydrogen Evolution Reaction.  Nano Lett. 2016;16:1097–1103.
  * [6] Wu L, Hofmann JP.  Comparing the Intrinsic HER Activity of Transition Metal Dichalcogenides: Pitfalls and Suggestions.  ACS Energy Lett. 2021;6:2619–2625.
  * [7] Yun Q, Li L, Hu Z, Lu Q, Chen B, Zhang H.  Layered Transition Metal Dichalcogenide-Based Nanomaterials for Electrochemical Energy Storage.  Advanced Materials. 2020;32(1903826).
  * [8] Tang B, Che B, Xu M, Ang ZP, Di J, Gao HJ, et al.  Recent Advances in Synthesis and Study of 2D Twisted Transition Metal Dichalcogenide Bilayers.  Small Structures. 2021;2(2000153).
  * [9] Cao Y, Robson AJ, Alharbi A, Roberts J, Woodhead CS, Noori YJ, et al.  Optical identification using imperfections in 2D materials.  2D Materials. 2017;4(045021).
  * [10] Wu Z, Luo Z, Shen Y, Zhao W, Wang W, Nan H, et al.  Defects as a factor limiting carrier mobility in WSe2: A spectroscopic investigation.  Nano Research. 2016;9:3622–3631.
  * [11] Chow PK, Jacobs-Gedrim RB, Gao J, Lu TM, Yu B, Terrones H, et al.  Defect-Induced Photoluminescence in Monolayer Semiconducting Transition Metal Dichalcogenides.  ACS Nano. 2015;9:1520–1527.
  * [12] Sumaiya SA, Liu J, Baykara MZ.  True Atomic-Resolution Surface Imaging and Manipulation under Ambient Conditions via Conductive Atomic Force Microscopy.  ACS Nano. 2022;16.
  * [13] Nowakowski K, Zandvliet HJW, Bampoulis P.  Barrier Inhomogeneities in Atomic Contacts on WS2.  Nano Lett. 2018;19:1190–1196.
  * [14] Severin N, Dzhanoev AR, Lin H, Rauf A, Kirstein S, Palma CA, et al.  Atomic resolution with high-eigenmode tapping mode atomic force microscopy.  Physical Review Research. 2022;4(023149).
  * [15] Barja S, Refaely-Abramson S, Schuler B, Qiu DY, Pulkin A, Wickenburg S, et al.  Identifying substitutional oxygen as a prolific point defect in monolayer transition metal dichalcogenides.  Nature Communications. 2019;10:1–8.
  * [16] Chang YS, Chen CY, Ho CJ, Cheng HR C M Chen, Fu TY, Huang YT, et al.  Surface electron accumulation and enhanced hydrogen evolution reaction in MoSe2 basal planes.  Nano Energy. 2021;84(105922).
  * [17] Addou R, Colombo L, Wallace RM.  Surface Defects on Natural MoS2.  ACS Appl Mater Interfaces. 2015;7:11921–11929.
  * [18] Blades WH, Frady NJ, Litwin PM, McDonnell SJ, Reinke P.  Thermally Induced Defects on WSe2.  J Phys Chem C. 2020;124:15337–15346.
  * [19] Yankowitz M, McKenzie D, LeRoy BJ.  Local Spectroscopic Characterization of Spin and Layer Polarization in WSe2.  Physical Review Letters. 2015;115(136803).
  * [20] Vancsó P, Magda GZ, Petö J, Noh JY, Kim YS, Hwang C, et al.  The intrinsic defect structure of exfoliated MoS2 single layers revealed by Scanning Tunneling Microscopy.  Scientific Reports. 2016;6(29726).
  * [21] Zhang C, Wang C, Yang F, Huang JK, Li LJ, Yao W, et al.  Engineering Point-Defect States in Monolayer WSe2.  ACS Nano. 2019;13:1595–1602.
  * [22] Mitterreiter E, Schuler B, Cochrane KA, Wurstbauer U, Weber-Bargioni A, Kastl C, et al.  Atomistic Positioning of Defects in Helium Ion Treated Single-Layer MoS2.  Nano Lett. 2020;20:4437–4444.
  * [23] Petö J, Ollár T, Vancsó P, Popov ZI, Magda GZ, Dobrik G, et al.  Spontaneous doping of the basal plane of MoS2 single layers through oxygen substitution under ambient conditions.  Nature Chemistry. 2018;10:1246–1251.
  * [24] Venugopal A, Pirkle A, Wallace RM, Colombo L, Vogel EM.  Contact Resistance Studies of Metal on HOPG and Graphene Stacks.  AIP Conference Proceedings. 2009;1173(324).
  * [25] Rosenberger MR, Chuang H, McCreary KM, Hanbicki AT, Sivaram SV, Jonker BT.  Nano-“Squeegee” for the Creation of Clean 2D Material Interfaces.  Appl Mater Interfaces. 2018;10:10379–10387.
  * [26] Bampoulis P, Van Bremen R, Yao Q, Poelsema B, Zandvliet HJW, Sotthewes K.  Defect Dominated Charge Transport and Fermi Level Pinning in MoS2/Metal Contacts.  ACS Applied Materials and Interfaces. 2017;9:19278–19286.
  * [27] Liu H, Han N, Zhao J.  Atomistic insight into the oxidation of monolayer transition metal dichalcogenides: from structures to electronic properties.  RSC Adv. 2015;5:17572-81.
  * [28] Longo RC, Addou R, KC S, Noh JY, Smyth CM, Barrera D, et al.  Intrinsic air stability mechanisms of two-dimensional transition metal dichalcogenide surfaces: basal versus edge oxidation.  2D Materials. 2017;4(2).
  * [29] Azcatl A, Qin X, Prakash A, Zhang C, Cheng L, Wang Q, et al.  Covalent Nitrogen Doping and Compressive Strain in MoS2 by Remote N2 Plasma Exposure.  Nano Letters. 2016;16:5437-43.
  * [30] González CB Biel, Dappe YJ.  Theoretical characterisation of point defects on a MoS2 monolayer by scanning tunnelling microscopy.  Nanotechnology. 2016;27(10):105702.
  * [31] Schuler B, Lee JH, Kastl C, Cochrane KA, Chen CT, Refaely-Abramson S, et al.  How Substitutional Point Defects in Two-Dimensional WS2 Induce Charge Localization, Spin–Orbit Splitting, and Strain.  ACS Nano. 2019;13:10520–10534.
  * [32] Zheng YJ, Chen Y, Huang YL, Gogoi PK, Li MY, Li LJ, et al.  Point Defects and Localized Excitons in 2D WSe2.  ACS Nano. 2019;13:6050–6059.
  * [33] Zhang S, Wang CG, Li MY, Huang D, Li LJ, Ji W, et al.  Defect Structure of Localized Excitons in a WSe2 Monolayer.  Phys Rev Lett. 2017;119(046101).
  * [34] Santosh K, Longo RC, Wallace RM, Cho K.  Surface oxidation energetics and kinetics on MoS2 monolayer.  Journal Of Applied Physics. 2015;117(135301).
  * [35] Iordanidou K, Houssa M, Pourtois G, Afanas’ev VV, A S.  Oxygen and hydroxyl adsorption on MS2 (M = Mo, W, Hf) monolayers: a first-principles molecular dynamics study.  Phys Status Solidi—RRL. 2016;10:787-91.
  * [36] Cochrane KA, Zhang T, Kozhakhmetov A, Lee JH, Zhang F, Dong C, et al.  Intentional carbon doping reveals CH as an abundant charged impurity in nominally undoped synthetic WS2 and WSe2.  2D Mater. 2020;7(3):031003.
  * [37] Cheng YC, Zhu ZY, Mi WB, Guo ZB, U S.  Prediction of two-dimensional diluted magnetic semiconductors: Doped monolayer MoS2 systems.  Phys Rev B. 2013;87(100401(R)).
  * [38] Rosenberger MR, Chuang HJ, McCreary KM, Li CH, Jonker BT.  Electrical Characterization of Discrete Defects and Impact of Defect Density on Photoluminescence in Monolayer WS2.  ACS Nano. 2018;12:1793–1800.
  * [39] Aghajanian M, Schuler B, Cochrane KA, Lee JH, Kastl C, Neaton JB, et al.  Resonant and bound states of charged defects in two-dimensional semiconductors.  Phys Rev B. 2020;101:081201(R).
  * [40] Mirabelli G, McGeough C, Schmidt M, McCarthy EK, Monaghan S, Povey IM, et al.  Air sensitivity of MoS2, MoSe2, MoTe2, HfS2, and HfSe2.  Journal of Applied Physics. 2016;120(125102).
  * [41] Ye F, Lee J, Hu J, Mao Z, Wei J, Feng PXL.  Environmental Instability and Degradation of Single- and Few-Layer WTe2 Nanosheets in Ambient Conditions.  Small. 2016;12:5802-8.
  * [42] Su TH, Lin YJ.  Effects of nitrogen plasma treatment on the electrical property and band structure of few-layer MoS2.  Applied Physics Letters. 2016;108(033103).
  * [43] Komsa HP, Kotakoski J, Kurasch S, Lehtinen O, Kaiser U, V KA.  Two-Dimensional Transition Metal Dichalcogenides under Electron Irradiation: Defect Production and Doping.  Physical Review Letters. 2012;109(035503).
  * [44] Cochrane KA, Zhang T, Kozhakhmetov A, Lee JH, Zhang F, Dong C, et al.  Intentional carbon doping reveals CH as an abundant charged impurity in nominally undoped synthetic WS2 and WSe2.  2D Materials. 2020;7(031003).
  * [45] Yurtsever A, ández Torre DF, González C, Jelínek P, Pou P, Sugimoto Y, et al.  Understanding image contrast formation in TiO2 with force spectroscopy.  Phys Rev B. 2012;85(125416).
  * [46] Wastl DS, Weymouth J, Giessibl FJ.  Optimizing atomic resolution of force microscopy in ambient conditions.  Phys Rev B. 2013;87(245415).
  * [47] Wastl DS, Michael J, Weymouth J, Giessibl FJ.  Atomic Resolution of Calcium and Oxygen Sublattices of Calcite in Ambient Conditions by Atomic Force Microscopy Using qPlus Sensors with Sapphire Tips.  ACS Nano. 2015;9:3858–3865.
  * [48] Schimmel T, Koch T, J K, M LS.  True atomic resolution under ambient conditions obtained by atomic force microscopy in the contact mode.  Applied Physics A. 1999;68:399–402.
  * [49] Korolkov VV, Summerfield A, Murphy A, Amabilino DB, Watanabe K, Taniguchi T, et al.  Ultra-high resolution imaging of thin films and single strands of polythiophene using atomic force microscopy.  Nature Communications. 2019;10(1537).
  * [50] Korolkov VV, Timokhin IG, Haubrichs R, Smith EF, Yang L, Yang S, et al.  Supramolecular networks stabilise and functionalise black phosphorus.  Nature Communications. 2017;8(1385).
  * [51] Bradford J, Cheng TS, James TSS, Khlobystov AN, Mellor CJ, Watanabe K, et al.  Graphene nanoribbons with hBN passivated edges grown by high-temperature molecular beam epitaxy.  2D Materials. 2023;10(3):035035.
  * [52] Enachescu M, Schleef D, Ogletree DF, Salmeron M.  Integration of point-contact microscopy and atomic-force microscopy: Application to characterization of graphite/Pt(111).  Phys Rev B. 1999;60(24).
  * [53] Celano U, Hantschel T, Giammaria G, Chintala RC, Conard T, Bender H, et al.  Evaluation of the electrical contact area in contact-mode scanning probe microscopy.  Journal of Applied Physics. 2015;117(214305).
  * [54] Custance O, Perez R, Morita S.  Atomic force microscopy as a tool for atom manipulation.  Nature Nanotechnology. 2009;4:pages 803–810.
  * [55] Pou P, Ghasemi SA, Jelinek P, Lenosky T, Goedecker S, Perez R.  Structure and stability of semiconductor tip apexes for atomic force microscopy.  Nanotechnology. 2009;20(264015).
  * [56] Jarvis SP, Kantorovich L, Moriarty P.  Structural development and energy dissipation in simulated silicon apices.  Beilstein J Nanotechnol. 2013;4:941–948.
  * [57] Sweetman A, Stirling J, Jarvis SP, Rahe P, Moriarty P.  Measuring the reactivity of a silicon-terminated probe.  Physical Review B. 2016;94(115440).
  * [58] Giessibl FJ.  Theory for an electrostatic imaging mechanism allowing atomic resolution of ionic crystals by atomic force microscopy.  Physical Review B. 1992;45(23).