# The cavity of 3CR 196.1: H$\alpha$ emission spatially associated with an
X-ray cavity

A. Jimenez-Gallardo Dipartimento di Fisica, Università degli Studi di Torino,
via Pietro Giuria 1, I-10125 Torino, Italy  European Southern Observatory,
Alonso de Córdova 3107, Vitacura, Región Metropolitana, Chile  INAF-
Osservatorio Astrofisico di Torino, via Osservatorio 20, 10025 Pino Torinese,
Italy  Dipartimento di Fisica e Astronomia dell’Università degli Studi di
Bologna, via P. Gobetti 93/2, I-40129 Bologna, Italy  E. Sani European
Southern Observatory, Alonso de Córdova 3107, Vitacura, Región Metropolitana,
Chile  F. Ricci Dipartimento di Fisica e Astronomia dell’Università degli
Studi di Bologna, via P. Gobetti 93/2, I-40129 Bologna, Italy  INAF-
Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, via Gobetti
93/3, I-40129 Bologna, Italy  C. Mazzucchelli Núcleo de Astronomía, Facultad
de Ingeniería y Ciencias, Universidad Diego Portales, Av. Ejército 441,
Santiago, 8320000, Chile  European Southern Observatory, Alonso de Córdova
3107, Vitacura, Región Metropolitana, Chile  B. Balmaverde INAF-Osservatorio
Astrofisico di Torino, via Osservatorio 20, 10025 Pino Torinese, Italy  F.
Massaro Dipartimento di Fisica, Università degli Studi di Torino, via Pietro
Giuria 1, I-10125 Torino, Italy  Istituto Nazionale di Fisica Nucleare,
Sezione di Torino, I-10125 Torino, Italy  INAF-Osservatorio Astrofisico di
Torino, via Osservatorio 20, 10025 Pino Torinese, Italy  Consorzio
Interuniversitario per la Fisica Spaziale, via Pietro Giuria 1, I-10125
Torino, Italy  A. Capetti INAF-Osservatorio Astrofisico di Torino, via
Osservatorio 20, 10025 Pino Torinese, Italy  W. R. Forman Center for
Astrophysics $|$ Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138,
USA  R. P. Kraft Center for Astrophysics $|$ Harvard & Smithsonian, 60 Garden
Street, Cambridge, MA 02138, USA  G. Venturi Instituto de Astrofísica,
Facultad de Física, Pontificia Universidad Católica de Chile, Casilla 306,
Santiago 22, Chile  INAF - Osservatorio Astrofisico di Arcetri, Largo E. Fermi
5, I-50125 Firenze, Italy  M. Gendron-Marsolais European Southern Observatory,
Alonso de Córdova 3107, Vitacura, Región Metropolitana, Chile  Instituto de
Astrofísica de Andalucía (IAA-CSIC), Glorieta de la Astronomía, 18008 Granada,
Spain  M. A. Prieto Departamento de Astrofísica, Universidad de La Laguna,
E-38206 La Laguna, Tenerife, Spain  Instituto de Astrofísica de Canarias
(IAC), E-38200 La Laguna, Tenerife, Spain  A. Marconi INAF - Osservatorio
Astrofisico di Arcetri, Largo E. Fermi 5, 50125 Firenze, Italy  Dipartimento
di Fisica e Astronomia, Università di Firenze, via G. Sansone 1, 50019 Sesto
F.no, Firenze, Italy  H. A. Peña-Herazo East Asian Observatory, 660 North
A’ohōkū Place, Hilo, HI 96720, USA  S. A. Baum  University of Manitoba, Dept.
of Physics and Astronomy, Winnipeg, MB R3T 2N2, Canada  C. P. O’Dea University
of Manitoba, Dept. of Physics and Astronomy, Winnipeg, MB R3T 2N2, Canada  L.
Lovisari INAF- Osservatorio di Astrofisica e Scienza dello Spazio di Bologna,
via Gobetti 93/3, I-40129 Bologna, Italy  R. Gilli INAF- Osservatorio di
Astrofisica e Scienza dello Spazio di Bologna, via Gobetti 93/3, I-40129
Bologna, Italy  E. Torresi INAF- Osservatorio di Astrofisica e Scienza dello
Spazio di Bologna, via Gobetti 93/3, I-40129 Bologna, Italy  A. Paggi
Dipartimento di Fisica, Università degli Studi di Torino, via Pietro Giuria 1,
I-10125 Torino, Italy  Istituto Nazionale di Fisica Nucleare, Sezione di
Torino, I-10125 Torino, Italy  INAF-Osservatorio Astrofisico di Torino, via
Osservatorio 20, 10025 Pino Torinese, Italy  V. Missaglia Dipartimento di
Fisica, Università degli Studi di Torino, via Pietro Giuria 1, I-10125 Torino,
Italy  Istituto Nazionale di Fisica Nucleare, Sezione di Torino, I-10125
Torino, Italy  INAF-Osservatorio Astrofisico di Torino, via Osservatorio 20,
10025 Pino Torinese, Italy  G. R. Tremblay Center for Astrophysics $|$ Harvard
& Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA  B. J. Wilkes Center
for Astrophysics $|$ Harvard & Smithsonian, 60 Garden Street, Cambridge, MA
02138, USA

###### Abstract

We present a multifrequency analysis of the radio galaxy 3CR 196.1
($z=0.198$), associated with the brightest galaxy of the cool core cluster
CIZAJ0815.4-0303. This nearby radio galaxy shows a hybrid radio morphology and
an X-ray cavity, all signatures of a turbulent past activity, potentially due
to merger events and AGN outbursts. We present results of the comparison
between Chandra and VLT/MUSE data for the inner region of the galaxy cluster,
on a scale of tens of kpc. We discovered H$\alpha$ \+ [N II]$\lambda 6584$
emission spatially associated with the X-ray cavity (at $\sim$10 kpc from the
galaxy nucleus) instead of with its rim. This result differs from previous
discoveries of ionized gas surrounding X-ray cavities in other radio galaxies
harbored in galaxy clusters and could represent the first reported case of
ionized gas filling an X-ray cavity, either due to different AGN outbursts or
to the cooling of warm ($10^{4}<T\leq 10^{7}$ K) AGN outflows. We also found
that the H$\alpha$, [N II]$\lambda\lambda 6548,6584$ and [S II]$\lambda\lambda
6718,6733$ emission lines show an additional redward component, at $\sim$1000
km s-1 from rest frame, with no detection in H$\beta$ or [O
III]$\lambda\lambda 4960,5008$. We believe the most likely explanation for
this redward component is the presence of a background gas cloud since there
appears to be a discrete difference in velocities between this component and
the rest frame.

galaxies: active — X-rays: general — radio continuum: galaxies — galaxies:
clusters: intracluster medium

##  1 Introduction

The extragalactic fraction of the Third Cambridge Catalog of radio sources and
its revised versions (3C, 3CR, 3CRR, Edge et al. 1959; Bennett 1962; Spinrad
et al. 1985; Laing et al. 1983) constitute one of the most valuable samples of
powerful radio-loud active galactic nuclei (AGN) to explore feedback processes
occurring between radio galaxies and their environments (see, e.g., McNamara &
Nulsen 2007, 2012, Fabian 2012 and Kraft et al. 2012).

One of the main pieces of evidence of radio mode feedback occurring in AGN in
galaxy clusters and groups is the presence of X-ray cavities, first reported
by Böhringer et al. (1996) and Churazov et al. (2000), due to the interaction
of jets with the intracluster medium (ICM; see reviews by McNamara & Nulsen
2007, Fabian 2012 and Gitti et al. 2012, and works by Bîrzan et al. 2004,
Cavagnolo et al. 2010, and McNamara & Nulsen 2012). Some of the most
remarkable examples of cavities include those of Perseus (see Fabian et al.
2003, 2006 and Graham et al. 2008), M87 (see e.g., Forman et al. 2005, 2017),
Hydra A (Nulsen et al. 2005a), Hercules A (Nulsen et al. 2005b), MS0735.6+7421
(McNamara et al. 2005, 2009), NGC 5813 (Randall et al. 2011, 2015), A2052
(Blanton et al. 2011), A2597 (Tremblay et al. 2012), and NGC 4636 (Jones et
al. 2002).

Although it is common to find H$\alpha$ emission either surrounding X-ray
cavities (see e.g., Blanton et al. 2011, Tremblay et al. 2015, Balmaverde et
al. 2018), or spatially associated with radio emission (known as the
“alignment effect”; see e.g., Fosbury 1986, Hansen et al. 1987, Baum et al.
1988, McCarthy 1988, Baum et al. 1990, Tremblay et al. 2009 and Baldi et al.
2019), the effect of radio jets on extended emission-line regions (EELR,
regions of line-emitting gas on scales of tens of kpc, see e.g., Baum et al.
1988) is not yet fully understood.

A recent analysis of radio, optical, and X-ray emission in 3CR 196.1 showed it
is a promising source to investigate how radio jets affect EELRs, thanks to
the discovery of an X-ray cavity at $\sim$10 kpc from the nucleus toward the
southwest (see Ricci et al. 2018). Moreover, recent observations from the
Multi-Unit Spectroscopic Explorer (MUSE; Bacon et al. 2010) at the Very Large
Telescope (VLT) obtained as part of the MUse RAdio Loud Emission line Snapshot
(MURALES) survey (Balmaverde et al. 2021) revealed that the H$\alpha$ \+ [N
II]$\lambda 6584$ emission extends beyond the radio emission and appears, in
projection, spatially associated with the X-ray cavity.

To test where such optical emission stands with respect to the X-rays, we
compared the optical, X-ray, and radio emission of 3CR 196.1 up to a few tens
of kpc, taking advantage of new VLT/MUSE observations. This emission had not
been detected in previous HST data originally published in Tremblay et al.
(2009), due to HST’s lower sensitivity compared to MUSE. This spatial
association of ionized gas with an X-ray cavity, instead of with its borders
has been rarely reported in the literature (see e.g., A1668, Pasini et al.
2021). We focus on the comparison of X-ray and optical observations,
previously investigated only separately by Ricci et al. (2018) and Balmaverde
et al. (2021), respectively.

This manuscript is organized as follows. A brief description of 3CR 196.1 is
given in § 2. Astrometric registration of all images is reported in § 3,
together with details about data reduction procedures for optical and X-ray
data sets. Results and their discussion are presented in § 4 and § 5,
respectively, while § 6 is devoted to our conclusions. Further details on the
astrometric registration and the ionized gas kinematics can be found in
Appendices A and B.

We adopted cgs units for numerical results and assumed a flat cosmology with
$H_{0}=69.6$ km s-1 Mpc-1, $\Omega_{M}=0.286$ and $\Omega_{\Lambda}=0.714$
(Bennett et al., 2014), unless otherwise stated. At the source redshift (i.e.,
$z=0.198$) the physical scale is 3.299 kpc/arcsec. Standard astronomical
orientation (north upwards and east to the left) is adopted throughout the
paper.

##  2 The strange case of 3CR 196.1

3CR 196.1 is the radio galaxy associated with the brightest cluster galaxy of
the galaxy cluster CIZA J0815.4-0303 (Kocevski et al. 2007). From an optical
perspective, this radio galaxy is classified as a low-excitation radio galaxy
(Buttiglione et al. 2010), while at radio frequencies it presents a hymor
radio structure (see Gopal-Krishna & Wiita 2000), with a classical edge-
darkened morphology, typical of FR Is, toward the southwest and a radio lobe,
typical of FR IIs on the opposite side (Fanaroff & Riley 1974).

Near-infrared (Madrid et al. 2006) and optical continuum observations (Baum et
al. 1988; de Koff et al. 1996) revealed that the host galaxy of 3CR 196.1
shows an elliptical morphology elongated from northeast to southwest, in the
same direction as the radio jet. Optical images also revealed the presence of
periodic shells possibly due to past merger events (Zirbel 1996), in agreement
with centroid shifts measured using optical isophotes of its host galaxy (de
Koff et al. 1996). Additionally, using HST emission line images, Baldi et al.
(2019) measured a lower limit for the total ionized gas mass for 3CR 196.1 of
3.5$\cdot$105 M⊙, using the nuclear electron density. The analysis of newly
obtained MUSE data focused on the [O III]$\lambda 5007$ properties by Speranza
et al. (2021) revealed a broad blue [O III]$\lambda 5007$ component with a
maximum velocity of $\sim$-640 km s-1, extending toward the northeast, in the
same direction as the radio jet.

Ricci et al. (2018) confirmed the presence of an X-ray surface brightness
depression, previously suggested by Massaro et al. (2012), by carrying out a
multiwavelength analysis of archival radio, optical, and X-ray observations.
This surface brightness depression indicates the presence of a “butterfly-
shaped” cavity at $\sim$10 kpc from the nucleus (see left panel of Fig. 1).
The temperature of the galaxy cluster ($\sim$4.2 keV) drops to $\sim$2.8 keV
in the inner 30 kpc, suggesting that 3CR 196.1 is the brightest cluster galaxy
of a cool core cluster.

At scales of tens of kpc, the H$\alpha$ \+ [N II]$\lambda 6584$ emission
observed with the Hubble Space Telescope is spatially aligned with the radio
jet and lies in the base of the X-ray cavity. Lastly, Ricci et al. concluded
that 3CR 196.1 could be an example where cold gas is being uplifted by radio
lobes and that it could have undergone several AGN outbursts at multiple
epochs. They also suggested its galaxy cluster may have experienced a merger
episode, revealed by the spiral pattern seen in the X-ray observation,
characteristic of gas sloshing (see, e.g., Markevitch & Vikhlinin 2007).

##  3 Data analysis procedures

X-ray data employed for this work are $\sim$8 ks archival Chandra observations
(ObsID 12729). The MUSE dataset was obtained from the MURALES survey (ID
0102.B-0048(A); two 10 minute exposures). Details on these observations as
well as on the Chandra and MUSE data reduction and analyses can be found in
Ricci et al. (2018) and Balmaverde et al. (2021), respectively.

###  3.1 Astrometric registration

Astrometric registration was carried out by comparing 8.4 GHz Very Large Array
(VLA111https://science.nrao.edu/facilities/vla/archive/index), Pan-
STARRS222https://catalogs.mast.stsci.edu/panstarrs/ $r$-band, MUSE, and
Chandra images. An image showing the sources used for the registration can be
found in Appendix A.

We decided against registering the Chandra image since the nucleus of 3CR
196.1 is not detected in the Chandra image and the typical shift applied to
Chandra images is $<$0.7″ (Massaro et al. 2011), which corresponds to a shift
of less than 2 pixels in our X-ray images and, thus, it is sufficiently
accurate for our comparison. Additionally, no registration was performed on
the radio since the radio core is not detected in the available radio maps and
their positional uncertainty is $<0.03\arcsec$ (adopting a conservative VLA
uncertainty of 10% of the beam
size333https://science.nrao.edu/facilities/vla/docs/manuals/oss/performance/positional-
accuracy). Thus, the final astrometric registration of MUSE data relied on the
Pan-STARRS $r$-band image since all sources in the MUSE data were also
detected in Pan-STARRS. A detailed description of the strategy adopted to
register the data can be found in Appendix A.

Figure 1: Left: Exposure-corrected 0.7 – 2 keV Chandra image of 3CR 196.1 with
H$\alpha$ \+ [N II]$\lambda 6584$ (blue) and 8.4 GHz VLA (cyan) contours. The
position of the centroid of the optical host, which corresponds to the radio
core position is shown as a black cross. The Chandra image has a pixel scale
of 0.492″ (corresponding to $\sim$1.6 kpc) and was smoothed with a 1.5″ radius
Gaussian kernel. X-ray contours are drawn at 0.25$\,\cdot\,$10-15 erg
cm${}^{-2}\,$s-1. H$\alpha$ \+ [N II]$\lambda 6584$ contours were drawn at 5,
10, 20, 50, and 100 times the rms. 8.4 GHz VLA contours were drawn at 5, 20,
and 50 times the rms level of the background. H$\alpha$ and the [N II]$\lambda
6584$ emission presents an offset with respect to the borders of the X-ray
cavity. The northeastern radio lobe is co-spatial with the X-ray emission
peak, while the southwestern lobe is co-spatial with the X-ray cavity. Right:
Redward [N II]$\lambda 6584$ emission with the same exposure-corrected 0.7 – 2
keV Chandra contours (red) and H$\alpha$ \+ [N II]$\lambda 6584$ contours
(blue) overlaid as in the left panel. The (a) and (b) regions marked in the
figure show the extraction areas of the spectra in Fig. 2

.

###  3.2 Chandra data analysis

We performed the X-ray data reduction following the standard procedures of the
Chandra Interactive Analysis of Observations (CIAO; Fruscione et al. 2006)
v4.11 threads444http://cxc.harvard.edu/ciao/threads/, adopting the Chandra
Calibration Database v4.8.4.1 (for an in-depth analysis of the X-ray
observation see Ricci et al. 2018).

X-ray flux maps were created by taking into account the exposure time and the
effective area. X-ray maps presented in this work were restricted to 0.7 - 2
keV to highlight the soft diffuse X-ray emission by using monochromatic
exposure maps set to the nominal energies of 1.2 keV. The flux map was
converted from units of counts s-1 cm-2 to cgs units by multiplying each event
by the nominal energy of the band, assuming that every event in the same band
has the same energy (see Massaro et al. 2015).

###  3.3 MUSE data analysis

Regarding the MUSE data analysis, we obtained the fully reduced data cube from
the ESO archive555http://archive.eso.org/scienceportal/home. We measured the
redshift by using the Penalized Pixel-Fitting code (Cappellari 2017) to fit
the stellar continuum and absorption features in a 2″ radius circular region
centered on the host position. Following this procedure, we measured a value
for the redshift of $z=0.1982$, which we adopted throughout our analysis.

####  3.3.1 Large scale ionized gas

We subtracted the continuum spaxel by spaxel from the archival MUSE reduced
data in two different ranges, blue range, from 5370 to 6170 $\AA$, and red
range, from 7170 to 8680 $\AA$, by fitting the continuum in each range using a
power law. We fitted the H$\beta$ and [O III]$\lambda\lambda 4960,5008$ lines
in the continuum subtracted blue range and the [O I]$\lambda 6300$, H$\alpha$,
[N II]$\lambda\lambda 6548,6584$ and [S II]$\lambda\lambda 6718,6733$ lines in
the continuum subtracted red range. In general, we fitted each line with a
single Gaussian component, fixing the wavelength separation between lines as
well as the line ratio of [N II]$\lambda\lambda 6548,6584$ and [O
III]$\lambda\lambda 4960,5008$ components, to their theoretical value of 1/3
(adopting CASE B, see Osterbrock & Ferland 2006).

####  3.3.2 Central region

Although a single Gaussian was adequate to fit the spectral components in most
spaxels, we found that the spectra in the central area and the X-ray cavity
region have additional components in H$\alpha$, [N II]$\lambda 6584$ and [S
II]$\lambda\lambda 6718,6733$. The position of the redward [N II]$\lambda
6584$ emission is shown in the right panel of Fig. 1.

To illustrate these redward components, we chose two regions in the central
area of 3CR 196.1, regions a and b, shown in the right panel of Fig. 1. Both
these regions were chosen to be close to the central region so the signal-to-
noise was the highest and with sizes large enough to match the seeing (i.e.,
1″), but not so large that the kinematics of the ionized gas would compromise
the identification of the different spectral components. Additionally, region
a was chosen to contain only a single spectral feature for each emission line,
while region b was chosen to illustrate a region with both “rest-frame” and
redward components.

A comparison of the spectra extracted from regions a and b is shown in Fig. 2,
where the rest-frame wavelengths of lines are marked in black and those of the
redward components are marked in red. Although there appear to be possible
additional spectral features corresponding to the H$\beta$ and [O I]$\lambda
6300$ lines, these components are not detected (with a significance below
2$\sigma$). The spectral position of the reward components in region b in Fig.
2 corresponds to $\sim$21 $\AA$ (i.e., $\sim$1000 km s-1) from rest-frame and
was derived on the basis of the spectral fit performed on the spectrum
extracted from region b. This fit, shown in Fig. 3, was performed using two
H$\alpha$ \+ [N II]$\lambda 6584$ triplets, one at rest-frame and one
redshifted by $\sim$21 $\AA$, and shows how the spectrum from region b cannot
be reproduced without introducing an additional redward component.

These redshifted spectral features are spectrally blended with emission at the
systemic redshift, i.e., the redward and “rest-frame” components have an
offset consistent with the separation between the H$\alpha$ and the [N
II]$\lambda 6584$ lines (see Fig. 2). This offset, together with the fact that
this redward component overlaps, in projection, with the AGN emission as well
as with the low signal-to-noise of the [N II]$\lambda 6548$ line, introduces a
degeneracy in the intensities of the H$\alpha$, [N II]$\lambda\lambda
6548,6584$ and [S II]$\lambda\lambda 6718,6733$ lines that we cannot break
without knowing a priori the ionization status of the source and the density
of the gas. Thus, the limited spectral ($\sim$50 km s-1) and spatial
($\sim$1″) resolutions of MUSE cubes prevent us from breaking the degeneracy,
and from providing a thorough analysis of the kinematics and ionization state
of the gas.

Figure 2: Spectra extracted from regions a and b (top and bottom panels in
Fig. 1, right panel, respectively). The top spectrum could be fitted using a
single Gaussian component for each line, while the bottom one shows an extra
redward component in all lines except for H$\beta$ and [O III]$\lambda\lambda
4960,5008$ lines. Black markers point to the rest-frame position of each line
while red markers are located at 21 $\AA$ from rest-frame ($\sim$1000 km s-1).
Uncertainties in the spectra are shown in gray. Figure 3: Fit to the spectrum
extracted from region b (i.e., a 1 arcsec2 region, as shown in the right panel
of Fig. 1) using a H$\alpha$ \+ [N II]$\lambda 6584$ at rest frame (yellow)
and another one redshifted by $\sim$21 $\AA$. Uncertainties in the spectrum
are shown in grey. Spectra in region b (see right panel of Fig. 1) cannot be
reproduced using a single triplet.

Figure 4: [O III]$\lambda\lambda 4960,5008$ (left) and [S II]$\lambda\lambda
6718,6733$ (right) emission with exposure-corrected 0.7 – 2 keV Chandra
contours (red) and H$\alpha$ \+ [N II]$\lambda 6584$ contours (blue) overlaid.
The position of the centroid of the optical host is indicated with a red
cross. The 1 arcsec2 extraction regions, L, C, F1, and F2, for the BPT diagram
in Fig. 5 are shown as black dashed squares. Chandra contours were smoothed
with a 1.5″ Gaussian kernel radius and drawn in red at 0.25$\,\cdot\,$10-15
erg cm${}^{-2}\,$s-1. H$\alpha$ \+ [N II]$\lambda 6584$ contours were drawn at
5, 10, 20, 50, and 100 times the rms. Region L marks the region co-spatial
with the northeastern radio lobe.

##  4 Results

We detected H$\alpha$+[N II]$\lambda 6584$ emission spatially associated with
the X-ray cavity (see left panel of Fig. 1) and, at the same time, did not
detect H$\beta$ or [O III]$\lambda 5007$ emission in the X-ray cavity region
(see left panel of Fig. 4). A comparison between X-ray and H$\alpha$+[N
II]$\lambda 6584$, [O III]$\lambda\lambda 4960,5008$, and [S
II]$\lambda\lambda 6718,6733$ is shown in the left panel of Fig. 1 and the
left and right panels of Fig. 4, respectively. Additionally, flux maps,
velocity, and velocity dispersion maps for all emission lines in the MUSE
observed range are shown in Fig. 8, Appendix B.

The morphology of the H$\alpha$+[N II]$\lambda 6584$ hints at differences in
the density of the environment. In particular, the northeastern component
(region L in the left panel of Fig. 4) being co-spatial with the northeastern
radio lobe and with the highest X-ray surface brightness peak indicates a
denser environment toward the northeast than toward the southwest, where the
X-ray cavity is located.

Figure 5: Location of 3CR 196.1 regions core, northeastern radio lobe and blue
and red filaments (C, L, F1 and F2, respectively) in the BPT diagnostic
diagrams. The size of regions was chosen to match the seeing of MUSE
observations (i.e., $\sim$1″). Extraction regions are shown in Fig. 4. The red
solid curves represent the Kewley et al. (2001) theoretical upper bound for
pure star formation, the red straight line shows the Kewley et al. (2006)
separation between AGN and LINERs, while the dashed red curve in the [N
II]$\lambda 6584$ BPT is the Kauffmann et al. (2003) empirical classification
separating star-forming galaxies and AGN. Contours represent the iso-densities
of all SDSS/DR7 emission line galaxies (Capetti & Baldi 2011). Extended
ionized gas structures in 3CR 196.1 have an ionization state compatible with
the LINER region of the BPT diagnostic diagrams, while the ionization state at
optical host position is compatible with Seyfert-like ionization.

We decided to give an overview of the gas kinematics via velocity channel
maps, as shown in Fig. 6 since the blending of the spectral components (see
Fig. 3) and the low signal-to-noise in the cavity region prevent us from
drawing strong quantitative conclusions from the study of the kinematics of
each component. This approach is commonly used for isolated lines. Thus, to
avoid contamination due to emission lines being part of a triplet, we
considered different reference lines for blue and red-shifted emission. We
considered as blue-shifted, emission blueward of the [N II]$\lambda 6548$ rest
frame line and, as red-shifted, emission redward of the [N II]$\lambda 6584$
rest frame. Therefore, the central wavelength range including the emission of
the H$\alpha$+[N II]$\lambda 6584$, dominated by spectral blending, is not
shown in the channel maps. Channel maps are shown in increments of 100 km s-1,
which correspond to two resolution elements, with the exception of the central
channels that are shown in increments of 50 km s-1.

Velocity channel maps show that the blue-shifted emission is co-spatial with
the north-eastern radio lobe (region L). Additionally, there seems to be a
filament of blue-shifted emission on the outer edge of the southern radio lobe
(F1). The red-shifted emission, up to $\sim$300 km s-1, shows some extended
emission that appears to follow the north-eastern radio jet and a filament
that seems to wrap around the southern lobe (F2). There also appears to be
some emission connecting filaments F1 and F2, which is most visible at
$\sim$400 km s-1 (A). This arched feature is apparently tracing the inner edge
of the X-ray cavity. Lastly, at the reddest velocities ($>$500 km s-1), the
only emission that is detected corresponds to the redward [N II]$\lambda 6584$
shown in the right panel of Fig. 1 and in the spectrum of region b in Fig. 2
(region R).

Figure 6: H$\alpha$ \+ [N II]$\lambda 6584$ flux in velocity bins of 100 km
s-1, with 8.4 GHz VLA (black) and exposure-corrected 0.7 - 2 keV Chandra
(blue) contours overlaid. MUSE images have a pixel size of 0.2 arcsec/pixel.
8.4 GHz VLA contours were drawn at 5, 20, and 50 times the rms level of the
background. Chandra contours were smoothed with a 1.5″ Gaussian kernel radius
and drawn in red at 0.25$\,\cdot\,$10-15 erg cm${}^{-2}\,$s-1. Centers of the
velocity bins are shown on the top right corners. Main features are labeled in
figures, with L, the region co-spatial with the northeastern radio lobe, F1
and F2, the blue and red filaments that appear to surround the southwestern
radio lobe, A, the arched emission that appears to trace the inner edge of the
X-ray cavity, F3, an additional filament not co-spatial with the X-ray cavity,
and R, the region polluted by the redward component.

All extended ionized gas features, i.e., lobe and red and blue-shifted
filaments (regions L, F1, and F2, marked in Fig. 4) present LINER-like
ionization states, as shown in the BPT diagram in Fig. 5 (see, e.g., Baldwin
et al. 1981 and Kewley et al. 2006). The only region shown in Fig. 5 with a
Seyfert-like ionization state is the region corresponding to the optical host
position. The position in the BPT diagram of the extra redward component shown
in Fig. 1, right panel, could not be established since the spectral and
spatial blending does not allow us to obtain estimates of the line ratios in
this region with the current signal-to-noise ratio of our MUSE data. All
extraction regions used for the BPT diagram were selected on the basis of
specific morphology/kinematics of the source and are shown in Fig. 4.
Additionally, all regions have a size of $\sim$1 arcsec2, to match the seeing
of the observation.

##  5 Discussion

Multiphase gas halos, surrounding AGN hosted in BCGs harbored in cool core
galaxy clusters, represent one of the keys to understanding AGN feedback
mechanisms. Currently, the main hypothesis for the origin of ionized gas in
these halos is the cooling and condensation of the hot, X-ray emitting plasma,
often in the shape of i) filaments (see e.g., Fabian et al. 2008, 2011; Fabian
2012, Gaspari et al. 2017, 2018, Qiu et al. 2020, 2021, Jimenez-Gallardo et
al. 2021), ii) rotating disks (see e.g., Wilman et al. 2005 Hamer et al.
2014), or iii) plumes (see e.g., Hamer et al. 2012, 2016, Pasini et al. 2019).
One of the most remarkable examples of ionized gas filaments in a cool core
cluster can be found in NGC 1275 in the center of the Perseus galaxy cluster
(Lynds 1970, Conselice et al. 2001 and Fabian et al. 2008), where ionized gas
filaments are spatially associated with an X-ray excess (Fabian et al. 2011).
The same situation occurs for other BCGs of cool core clusters, where a tight
connection between X-ray and optical filamentary emission was found, such as
the cases of A1795 (Crawford et al. 2005), A1644 (McDonald et al. 2010), A2597
(Tremblay et al. 201), and MKW 3s (Jimenez-Gallardo et al. 2021).

Ionized gas filaments have also been found surrounding X-ray cavities, i.e.,
co-spatial with the rims of X-ray cavities (see e.g., works on A2052 by
Blanton et al. 2011 and Balmaverde et al. 2018, as well as works by Lim et al.
2008, Tremblay et al. 2015 and Olivares et al. 2019). However, ionized gas
spatially associated, and thus, potentially filling an X-ray cavity, is a rare
occurrence with no clear explanation. In this work, we present the potential
first detection of ionized gas spatially associated with an X-ray cavity,
using H$\alpha$ \+ [N II]$\lambda 6584$ emission as a tracer. Although the
scenario where this ionized gas is actually overlaid to the X-ray cavity
cannot be conclusively excluded, it appears unlikely. If that was the case, we
would expect to see ionized gas spatially associated with the outer borders of
the X-ray cavity, where its column density would be the highest, instead of
with the X-ray cavity.

###  5.1 The north-eastern plume

Similar to other cool core galaxy clusters such as A1991, A3444, Ophiuchus
(Hamer et al. 2012), A111, A133, A2415 (Hamer et al. 2016), and A2495 (Pasini
et al. 2019), the peak of the X-ray emission surrounding 3CR 196.1 does not
correspond to the position of the BCG, presenting an offset of $\sim$2″ (i.e.,
$\sim$6.5 kpc; see left panel of Fig. 2). Such offsets have been interpreted
as caused by sloshing of the ICM due to events such as minor mergers.
Typically, in those cases, the ionized gas forms plumes either following the
X-ray peak or connecting the X-ray peak and the BCG. However, in the case of
3CR 196.1, we argue that the X-ray peak seen in the ICM is most likely caused
by the interaction of the north-eastern radio jet with the ICM, as seen in 3CR
171 and 3CR 305 (Hardcastle et al. 2010, 2012), due to the alignment between
the north-eastern radio lobe and the X-ray peak, as well as to the “hammer-
head” shape of the north-eastern radio lobe and the ionized gas plume
extending in the direction opposite to the position of the BCG (see Fig. 4).
Thus, the ionized gas corresponding to region L could originate from the
cooling of the jet-heated ICM. This origin is compatible with the ionized gas
presenting a LINER-like ionization state, as shown in Fig. 5 since LINER-like
emission in extranuclear regions is often associated with large-scale outflows
and shocks or with radio jet-shocked regions (see Kewley et al. 2006 and
references therein).

###  5.2 Ionized gas filaments around 3CR 196.1

One of the most likely mechanisms responsible for the formation of ionized gas
filaments in cool core galaxy clusters is chaotic cold accretion (CCA, see
Gaspari et al. 2012, 2013, 2015, 2017, 2018), i.e., the condensation of the
hot plasma in the ICM and its precipitation onto the black hole. ICM cooling
due to non-linear thermal instabilities was first introduced by Pizzolato &
Soker (2005) as “cold feedback”. More recently, Gaspari et al. (2018) adopted
the ratio of the cooling time and the eddy turnover time,
$C=t_{cool}/t_{eddy}$, as a tool to assess the multiphase state of a system,
with $0.6<C<1.8$ marking the extent of the condensation region. We applied
this criterion to 3CR 196.1, where we assumed $t_{cool}\sim$500 Myr, as
obtained by Ricci et al. (2018) for the inner region (i.e., 10″ or $\sim$30
kpc radius), since currently available X-ray observations did not allow us to
carry out a more detailed spectral analysis. We obtained the eddy turnover
time, as a function of the distance to the galaxy cluster core, as
$t_{eddy}=2\pi\frac{r^{2/3}L^{1/3}}{\sigma_{v,L}}$ (see Gaspari et al. 2018),
where $L\sim$10 kpc is the injection scale, traced by the diameter of the
X-ray cavity, and $\sigma_{v,L}\sim$212 km s-1 is the velocity dispersion at
the injection scale. Thus, significant precipitation is expected to occur
within a region of $r\sim$9 - 48 kpc radius, where $C\approx 1$. Since the
blue and red filaments at the edges of the southern radio lobe (F1 and F2 in
the right panel of Fig. 4) extend up to 10 kpc (where $C\sim$1.7) and the
small ionized gas filament, F3, extends up to 15 kpc (where $C\sim$1.3) to the
south-east (see Fig. 6), the origin of ionized gas filaments surrounding 3CR
196.1 is compatible with a condensation rain. Moreover, the velocity and
velocity dispersion along filament F2, $\sim$200 km s-1 and $\sim$100 km s-1,
respectively, are consistent with those literature observations compared with
CCA simulations by Gaspari et al. (2018). Therefore, 3CR 196.1 presents
possible signs of black hole feeding.

Alternatively, as proposed e.g. in simulations by Qiu et al. (2020, 2021),
ionized gas filaments in 3CR 196.1 could be due to AGN outflows uplifting warm
($10^{4}<T\leq 10^{7}$ K) gas from the central 2 kpc of the galaxy cluster,
which then experiences radiative cooling. The energy needed to inflate
adiabatically an ionized gas bubble with a radius $r\sim$5 kpc, i.e., the size
of the arched structure A in Fig. 6, and velocity $v\sim$400 km s-1, in an
ambient medium with density $n_{0}=0.098$ cm-3, as obtained by Ricci et al.
(2018) for the inner $\sim$30 kpc of 3CR 196.1, can be obtained following the
approach described by Nesvadba et al. (2006):

$\dot{E}\approx 1.5\cdot
10^{46}\,r^{2}_{10}\,v^{3}_{1000}\,n_{0}\,\text{ergs}^{-1}$ (1)

where $r_{10}$ is the radius of the bubble in units of 10 kpc, and $v_{1000}$
is the velocity in units of 1000 km s-1. Thus, the minimum energy needed to
inflate such a bubble, $\dot{E}\sim$2.4$\cdot 10^{43}$ erg s-1, is an order of
magnitude below the jet power obtained by Ricci et al. (2018),
$P_{jet}\sim$1.9$\cdot 10^{44}$ erg s-1, and, therefore, the ionized gas
spatially associated with the X-ray cavity is consistent with a bubble being
inflated by the radio jet.

The electron density in filaments F1 and F2 can be estimated using the [S
II]$\lambda\lambda 6718,6733$ duplet ratio as described in Sanders et al.
(2016). We selected two regions of the size of the PSF ($\sim$1.0″) over the
filaments to compute the electron density. Thus, the electron density within
the cavity is $<350$ e${}^{-}\,$cm-3. This density is comparable with that of
the filaments found surrounding an X-ray cavity in A2052 by Balmaverde et al.
(2018), which were interpreted as being caused by an expanding bubble, as well
as with other filaments in cool core galaxy clusters (see e.g., McDonald et
al. 2012, Jimenez-Gallardo et al. 2021). Additionally, extended ionized gas
features, corresponding to filaments F1 and F2, could be interpreted as being
caused by the interaction of the southern radio jet with its surrounding
medium, since the filaments present LINER-like ionization states, as shown in
Fig. 5.

Lastly, similarly to the expanding bubble in A2052, filaments F1 and F2 in 3CR
196.1 are blue and redshifted, respectively, and a potential region with split
lines that could be due to the overlap of the emission of both filaments, a
typical sign of an expanding bubble (Balmaverde et al. 2018), is also found to
be co-spatial with the X-ray cavity (region R in Fig. 6). If that is the case,
the expansion velocity of the bubble, which can only be estimated indirectly
through X-ray observations, could be measured as the maximum velocity observed
at the arched structure, A, i.e., $\sim$400 - 500 km s-1. Nevertheless, due to
the low signal-to-noise ratio in the cavity region of current MUSE
observations, the kinematics of the ionized gas filaments cannot be firmly
derived and, thus, this scenario cannot be confirmed.

###  5.3 The central region of 3CR 196.1

During this analysis, we discovered that, in the central and X-ray cavity
regions, all spectral features, except for the H$\beta$ and [O
III]$\lambda\lambda 4960,5008$ and [O I]$\lambda 6300$ lines, have an
additional component at $\sim$1000 km s-1 to the red of rest-frame (see right
panel of Fig. 1 and Figs. 2 and 3). This additional redward component could be
due to a background gas cloud since there does not appear to be a continuous
range of velocities from rest-frame to the redward component that would
indicate the presence of an outflow. This emission could also be due to
filaments of ionized gas being accreted onto the black hole through CCA,
similarly to filament F3. Despite that, the spatial resolution ($\sim$1″) of
the current observations together with the blending of lines in the H$\alpha$
\+ [N II]$\lambda 6584$ complex due to the limited spectral resolution
($\sim$50 km s-1) prevent an in-depth analysis on the kinematics and
ionization state of the gas in the central and X-ray cavity regions with
current signal-to-noise.

##  6 Summary and conclusions

We carried out a comparison of optical VLT/MUSE and X-ray Chandra observations
of 3CR 196.1, the central radio galaxy associated with the cool core galaxy
cluster CIZA J0815.4-0303, which exhibits a “butterfly-shaped” X-ray cavity at
$\sim$10 kpc toward the southwest of the nucleus. Through this analysis, we
detected the presence of ionized gas spatially associated with an X-ray
cavity.

The main possible explanations for detecting H$\alpha$ \+ [N II]$\lambda 6584$
emission spatially associated with an X-ray cavity instead of with its rim
include:

  1. 1.

The ionized gas is actually concentrated in filaments that wrap around the
X-ray cavity, such as the ones shown in Fig. 6. However, due to projection
effects, as well as to the presence of an additional redward component in the
cavity region (see right panel of Fig. 1 and Fig. 2), the ionized gas appears
to be co-spatial with the X-ray cavity. We tend to disfavor this scenario
since it implies that we would be seeing the gas where it has the lowest
column density along the line of sight, instead of at the edges, where its
column density would be the highest.

  2. 2.

The ionized gas is actually filling the X-ray cavity. In that case, there
could be two possible explanations: (i) the gas has undergone multiple
ionization events due to different AGN outbursts, or (ii) the ionized gas is
forming filaments originating from the cooling of warm ($10^{4}<T\leq 10^{7}$
K) AGN outflows. Scenario (i) is compatible with an AGN outburst occurring
after the one responsible for the X-ray cavity, since typical cycling times
between quiescence and activity in radio galaxies are $\sim$107-108 yr (see
Shabala et al. 2008 and references therein) and the estimated age of the X-ray
cavity is 12 Myr (Ricci et al. 2018). On the other hand, scenario (ii) is
compatible with simulations by Qiu et al. (2020), which show that warm
($10^{4}<T\leq 10^{7}$ K) gas present in the central 2 kpc of galaxy clusters
(due to radiative feedback) could be lifted by AGN outflows and radiatively
cooled, originating ionized gas filaments. Additional simulations in Qiu et
al. (2021) showed that these filaments could be shaped as rings perpendicular
to the direction of the outflow, which would explain the arched feature (A in
Fig. 6), that appears to connect filaments F1 and F2 (surrounding the
southwestern radio lobe in Fig. 4, right panel) and that seems to trace the
inner edge of the X-ray cavity.

To distinguish between these scenarios, additional observations would allow us
to characterize, in more detail, the temperature and density profiles of the
central ICM emission. To characterize the ionized gas, we would also need
deeper and higher spatial resolution MUSE observations, with, e.g., the
adaptive optic-assisted MUSE narrow-field-mode. Alternatively, the new ERIS
instrument at the VLT, or the NIRSPEC IFU on the James Webb Space Telescope,
could help us understand the origin of the additional redward component found
in the spectra, as well as to derive the kinematics of the ionized gas
spatially associated with the X-ray cavity. Such observations would allow us
to measure the expansion velocity in the case in which the ionized gas is
tracing an expanding bubble.

## References

  * Bacon et al. (2010) Bacon, R., Accardo, M., Adjali, L., et al. 2010, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 7735, Proc. SPIE, 773508
  * Baldi et al. (2019) Baldi, R. D., Rodríguez Zaurín, J., Chiaberge, M., et al. 2019, ApJ, 870, 53
  * Baldwin et al. (1981) Baldwin, J. A., Phillips, M. M. & Terlevich, R. 1981, PASP, 93, 5
  * Balmaverde et al. (2018) Balmaverde, B., Capetti, A., Marconi, A., et al. 2018, A&A, 612, A19
  * Balmaverde et al. (2019) Balmaverde, B., Capetti, A., Marconi, A., et al. 2019, A&A, 632, A124
  * Balmaverde et al. (2021) Balmaverde, B., Capetti, A., Marconi, A., et al. 2021, A&A, 645, A12
  * Baum et al. (1988) Baum, S. A., Heckman, T. M., Bridle, A., et al. 1988, ApJS, 68, 643
  * Baum et al. (1990) Baum, S. A., Heckman, T. & van Breugel, W. 1990, ApJS, 74, 389
  * Bennett (1962) Bennett, A. S. 1962, MNRAS, 125, 75
  * Bennett et al. (2014) Bennett, C. L., Larson, D., Weiland, J. L. & Hinshaw, G. 2014, ApJ, 794, 135
  * Bîrzan et al. (2004) Bîrzan, L., Rafferty, D. A., McNamara, B. R., et al. 2004, ApJ, 607, 800
  * Blanton et al. (2011) Blanton, E. L., Randall, S. W., Clarke, T. E., et al. 2011, ApJ, 737, 99
  * Böhringer et al. (1996) Böhringer, H., Neumann, D. M., Schindler, S., et al. 1996, ApJ, 467, 168
  * Buttiglione et al. (2010) Buttiglione, S., Capetti, A., Celotti, A., et al. 2010, A&A, 509, A6
  * Capetti & Baldi (2011) Capetti, A. & Baldi, R. D. 2011, A&A, 529, A126
  * Cappellari (2017) Cappellari, M. 2017, MNRAS, 466, 798
  * Cavagnolo et al. (2010) Cavagnolo, K. W., McNamara, B. R., Nulsen, P. E. J., et al. 2010, ApJ, 720, 1066
  * Churazov et al. (2000) Churazov, E., Forman, W., Jones, C., et al. 2000, A&A, 356, 788
  * Conselice et al. (2001) Conselice, C. J., Gallagher, J. S., III & Wyse, R. F. G. 2001, AJ, 122, 2281
  * Crawford et al. (2005) Crawford, C. S., Sanders, J. S. & Fabian, A. C. 2005, MNRAS, 361, 17
  * de Koff et al. (1996) de Koff, S., Baum, S. A., Sparks, W. B., et al. 1996, ApJS, 107, 621
  * Edge et al. (1959) Edge, D. O., Shakeshaft, J. R., McAdam, W. B., et al. 1959, MmRAS, 69, 37
  * Fabian et al. (2003) Fabian, A. C., Sanders, J. S., Allen, S. W., et al. 2003, MNRAS, 344, L43
  * Fabian et al. (2006) Fabian, A. C., Sanders, J. S., Taylor, G. B., et al. 2006, MNRAS, 366, 417
  * Fabian et al. (2008) Fabian, A. C., Johnstone, R. M., Sanders, J. S., et al. 2008, Nature, 454, 968
  * Fabian et al. (2011) Fabian, A. C., Sanders, J. S., Williams, R. J. R., et al. 2011, MNRAS, 417, 172
  * Fabian (2012) Fabian, A. C. 2012, ARA&A, 50, 455
  * Fanaroff & Riley (1974) Fanaroff, B. L. & Riley, J. M. 1974, MNRAS, 167, 31P
  * Forman et al. (2005) Forman, W., Jones, C., Churazov, E., et al. 2005, ApJ, 665, 1057
  * Forman et al. (2017) Forman, W., Churazov, E., Jones, C., et al. 2017, ApJ, 844, 122
  * Fruscione et al. (2006) Fruscione, A., McDowell, J. C., Allen, G. E., et al. 2006, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, 6270, 62701V
  * Fosbury (1986) Fosbury, R. A. E. 1986, in Astrophysics and Space Science Library Vol. 121, Structure and Evolution of Active Galactic Nuclei, ed. G. Giuricin, M. Mezzetti, M. Ramella, & F. Mardirossian (Dordrecht: Reidel), 297
  * Gaspari et al. (2012) Gaspari, M., Ruszkowski, M. &Sharma, P. 2012, ApJ, 746, 94
  * Gaspari et al. (2013) Gaspari, M., Ruszkowski, M. & Oh, S. Peng 2013, MNRAS, 432, 3401
  * Gaspari et al. (2015) Gaspari, M., Brighenti, F. & Temi, P. 2015, A&A, 579, A62
  * Gaspari et al. (2017) Gaspari, M., Temi, P. & Brighenti, F. 2017, MNRAS, 466, 677
  * Gaspari et al. (2018) Gaspari, M., McDonald, M., Hamer, S. L., et al. 2018, ApJ, 854, 167
  * Gitti et al. (2012) Gitti, M., Brighenti, F. & McNamara, B. R. 2012, AdAst, 2012, 950641
  * Gopal-Krishna & Wiita (2000) Gopal-Krishna & Wiita, P. J. 2000, A&A, 363, 507
  * Graham et al. (2008) Graham, J., Fabian, A. C. & Sanders, J. S. 2008, MNRAS, 386, 278
  * Hamer et al. (2012) Hamer, S. L., Edge, A. C., Swinbank, A. M., et al. 2012, MNRAS, 421, 3409
  * Hamer et al. (2014) Hamer, S. L., Edge, A. C., Swinbank, A. M., et al. 2014, MNRAS, 437, 862
  * Hamer et al. (2016) Hamer, S. L., Edge, A. C., Swinbank, A. M., et al. 2016, MNRAS, 460, 1758
  * Hansen et al. (1987) Hansen, L., Norgaard-Nielsen, H. U. & Jorgensen, H. E. 1987, A&AS, 71, 465
  * Hardcastle et al. (2010) Hardcastle, M. J., Massaro, F. & Harris, D. E., 2010, MNRAS, 401, 2697
  * Hardcastle et al. (2012) Hardcastle, M. J., Massaro, F. & Harris, D. E., 2012, MNRAS, 424, 1774
  * Jimenez-Gallardo et al. (2021) Jimenez-Gallardo, A., Massaro, F., Balmaverde, B., et al. 2021, ApJ, 912, L25
  * Jones et al. (2002) Jones, C., Forman, W.,Vikhlinin, A., et al. 2002, ApJ, 567, L115
  * Kauffmann et al. (2003) Kauffmann, G., Heckman, T. M., Tremonti, C., et al. 2003, MNRAS, 346, 1055
  * Kewley et al. (2001) Kewley, L. J., Heisler, C. A., Dopita, M. A, et al. 2001, ApJS, 132, 37
  * Kewley et al. (2006) Kewley, L. J., Groves, B., Kauffmann, G., et al. 2006, MNRAS, 372, 961
  * Kocevski et al. (2007) Kocevski, D. D., Ebeling, H., Mullis, C. R., et al. 2007, ApJ, 662, 224
  * Kraft et al. (2012) Kraft, R. P., Birkinshaw, M., Nulsen, P. E. J., et al. 2012, ApJ, 749, 19
  * Laing et al. (1983) Laing, R. A., Riley, J. M. & Longair, M. S. 1983, MNRAS, 204, 151
  * Lim et al. (2008) Lim, J., Ao, YiPing & Dinh-V-Trung 2008, ApJ, 672, 252
  * Lynds (1970) Lynds, R. 1970, ApJ, 159, L151
  * Madrid et al. (2006) Madrid, J. P., Chiaberge, M., Floyd, D., et al. 2006, ApJS, 164, 307
  * Markevitch & Vikhlinin (2007) Markevitch, M. & Vikhlinin, A. 2007, Phys. Rep., 443, 1
  * Massaro et al. (2011) Massaro, F., Harris, D. E. & Cheung, C. C. 2011, ApJS, 197, 24
  * Massaro et al. (2012) Massaro, F., Tremblay, G. R., Harris, D. E., et al. 2012, ApJS, 203, 31
  * Massaro et al. (2015) Massaro, F., Harris, D. E., Liuzzo, E., et al. 2015, ApJS, 220, 5
  * McCarthy (1988) McCarthy, P. J. 1988, PhD thesis, AA (California Univ., Berkeley)
  * McDonald et al. (2010) McDonald M., Veilleux, S., Rupke, D. S. N., et al. 2010, ApJ, 721, 1262
  * McDonald et al. (2012) McDonald M., Veilleux, S., Rupke, D. S. N., et al. 2012, ApJ, 746, 153
  * McNamara et al. (2005) McNamara, B. R., Nulsen, P. E. J., Wise, M. W., et al. 2005, Nature, 433, 7021
  * McNamara & Nulsen (2007) McNamara, B. R. & Nulsen, P. E. J. 2007, ARA&A, 45, 117
  * McNamara et al. (2009) McNamara, B. R., Kazemzadeh, F., Rafferty, D. A., et al. 2009, ApJ, 698, 594
  * McNamara & Nulsen (2012) McNamara B. R. & Nulsen, P. E. J. 2012, New Journal of Physics, 14, 5
  * Nesvadba et al. (2006) Nesvadba, N. P. H., Lehnert, M. D., Eisenhauer, F., et al. 2006, ApJ, 650, 693
  * Nulsen et al. (2005a) Nulsen, P. E. J., McNamara, B. R., Wise, M. W., et al. 2005a, ApJ, 628, 629
  * Nulsen et al. (2005b) Nulsen, P. E. J., Hambrick, D. C., McNamara, B. R., et al. 2005b, ApJ, 625, L9
  * Olivares et al. (2019) Olivares, V., Salome, P., Combes, F., et al. 2019, A&A, 631, A22
  * Osterbrock & Ferland (2006) Osterbrock, D. E. & Ferland, G. J. 2006, Astrophysics of gaseous nebulae and active galactic nuclei (2nd ed.; Sausalito, CA: Science Books)
  * Pasini et al. (2019) Pasini, T., Gitti, M., Brighenti, F., et al. 2019, ApJ, 885, 111
  * Pasini et al. (2021) Pasini, T., Gitti, M., Brighenti, F., et al. 2021, ApJ, 911, 66
  * Pizzolato & Soker (2005) Pizzolato, F. & Soker, N. 2005, ApJ, 632, 821
  * Prieto et al. (2016) Prieto, J. L., Krühler, T., Anderson, J. P., et al. 2016, ApJ, 830, L32
  * Qiu et al. (2020) Qiu, Y., Bogdanović, T., Li, Y., et al. 2020, Nature Astronomy, 4, 900
  * Qiu et al. (2021) Qiu, Y., Hu, H., Inayoshi, K., et al. 2021, ApJ, 917, L7
  * Randall et al. (2011) Randall, S. W., Forman, W. R., Giacintucci, S., et al. 2011, ApJ, 726, 86
  * Randall et al. (2015) Randall, S. W., Nulsen, P. E. J., Jones, C., et al. 2015, ApJ, 805, 112
  * Ricci et al. (2018) Ricci, F., Lovisari, L., Kraft, R. P., et al. 2018, ApJ, 867, 35
  * Sanders et al. (2016) Sanders, Ryan L., Shapley, A. E., Kriek, M., et al. 2016, ApJ, 825, L23
  * Shabala et al. (2008) Shabala, S. S., Ash, S., Alexander, P., et al. 2008, MNRAS, 388, 625
  * Speranza et al. (2021) Speranza, G., Balmaverde, B., Capetti, A., et al. 2021, A&A, 653, A150
  * Spinrad et al. (1985) Spinrad, H., Djorgovski, S., Marr, J. & Aguilar, L. 1985, PASP, 97, 932
  * Tremblay et al. (2009) Tremblay, G. R., Chiaberge, M., Sparks, W. B., et al. 2009, ApJS, 183, 278
  * Tremblay et al. (2012) Tremblay, G. R., O’Dea, C. P., Baum, S. A., et al. 2012, MNRAS, 424, 1026
  * Tremblay et al. (2015) Tremblay, G. R., O’Dea, C. P., Baum, S. A., et al. 2015, MNRAS, 451, 3768
  * Tremblay et al. (201) Tremblay, G. R., Oonk, J. B. R., Combes, F., et al. 2016, Nature, 534, 218
  * Tucker et al. (2021) Tucker, M. A., Shappee, B. J., Hinkle, J. T., et al. 2021, MNRAS, 506, 6014
  * Wilman et al. (2005) Wilman, R. J., Edge, A. C. & Johnstone, R. M. 2005, MNRAS, 359, 755
  * Zirbel (1996) Zirbel, E. L. 1996, ApJ, 473, 713

We thank the anonymous referee for their useful comments that led to the
substantial improvement of the paper. A. J. thanks M. Gaspari for his feedback
on Chaotic Cold Accretion. F. M. is in debt with S. Bianchi for his valuable
input on X-ray photoionization scenarios. This work is supported by the
“Departments of Excellence 2018 - 2022” Grant awarded by the Italian Ministry
of Education, University and Research (MIUR) (L. 232/2016). This research has
made use of resources provided by the Ministry of Education, Universities and
Research for the grant MASF_FFABR_17_01. This investigation is supported by
the National Aeronautics and Space Administration (NASA) grants GO9-20083X and
GO0-21110X. A. J. acknowledges the financial support (MASF_CONTR_FIN_18_01)
from the Italian National Institute of Astrophysics under the agreement with
the Instituto de Astrofisica de Canarias for the “Becas Internacionales para
Licenciados y/o Graduados Convocatoria de 2017”. F.R. acknowledges support
from PRIN MIUR 2017 project “Black Hole winds and the Baryon Life Cycle of
Galaxies: the stone-guest at the galaxy evolution supper”, contract
$\\#$2017PH3WAT. G.V. acknowledges support from ANID program FONDECYT
Postdoctorado 3200802. A.P. acknowledges financial support from the Consorzio
Interuniversitario per la fisica Spaziale (CIFS) under the agreement related
to the grant MASF_CONTR_FIN_18_02. W.F. and R.K. acknowledge support from the
Smithsonian Institution and the Chandra High Resolution Camera Project through
NASA contract NAS8-03060. S.B. and C.O. acknowledge support from the Natural
Sciences and Engineering Research Council (NSERC) of Canada.

##  Appendix A Atrometric registration

All details related to the strategy adopted to carry out the astrometric
registration of MUSE images are reported in the following. We measured the
centroids of the optical counterpart of 3CR 196.1 and sources SRC1 to SRC4, in
Fig. 7, in the Pan-STARRS and MUSE white light collapsed images and applied
the average of the shifts to the MUSE data cube. The final shift was 3.1″
(corresponding to $\sim$10 kpc, which is expected since it is a known issue
that MUSE astrometry can be off by several arcseconds, as reported by Prieto
et al. 2016, Balmaverde et al. 2019, 2021, and Tucker et al. 2021). We then
verified that the registration adopted was appropriate for sources SRC5 to
SRC8, in Fig. 7, by measuring the difference between their centroids in the
Pan-STARRS and MUSE white light collapsed images after the shift. The adopted
shift yielded an rms of 0.18″. The final alignment between Pan-STARRS and MUSE
data for 3CR 196.1 and the rest of the field sources is shown in Fig. 7, for
the MUSE white light image (red contours) and the H$\alpha$+[N II]$\lambda
6584$ emission (blue contours).

Figure 7: r-band PanSTARRS image with MUSE 4800 - 6800 $\AA$ (red) and MUSE
H$\alpha$ \+ [N II]$\lambda 6584$ (blue) contours overlaid. Sources SRC1 to
SRC4 were used to compute the shift applied, while sources SRC5 to SRC8
allowed us to verify the registration. MUSE images have a pixel size of 0.2″.
4800 - 6800 $\AA$ contours were drawn at 20, 50, and 200 times the root mean
square (rms) level of the background. H$\alpha$+[N II]$\lambda 6584$ contours
were drawn at 5, 10, 20, 50, and 100 times the rms level of the background. We
found a good alignment between the Pan-STARRS and the MUSE emission applying a
shift o 3.1″ (i.e., $\sim$10 kpc).

##  Appendix B Ionized gas kinematics

Here we report flux, velocity, and velocity dispersion maps for the H$\beta$,
[O III]$\lambda\lambda 4960,5008$, [O I]$\lambda 6300$, H$\alpha$, [N
II]$\lambda\lambda 6548,6584$, and [S II]$\lambda\lambda 6718,6733$ emission
lines (see Fig. 8), obtained following the fitting procedure described in §
3.3, with Chandra and 8.4 GHz contours overlaid in blue and black,
respectively. All maps were spatially binned by a factor of 2, resulting in a
pixel size of 0.4 arcsec/pixel, to spatially Nyquist-sample the observation,
as the image quality of MUSE data is of the order of 1″. Additionally this
binning allows us to better match the resolution of Chandra and to increase
the signal-to-noise ratio.

All emission features show similar morphologies and kinematics. The largest
blueshift ($\sim$-400 km s-1) can be seen at the termination of the
northeastern radio lobe, pointing to the cooling of jet-heated ICM as the
origin of ionized gas in the north-east direction. Lastly, the velocity
dispersion of all lines is higher ($\sim$200 - 300 km s-1) in the central
$\sim$10 kpc of 3CR 196.1.

Figure 8: Flux (left), velocity (middle), and velocity dispersion (right) maps
with 8.4 GHz VLA (black) and exposure-corrected 0.7 - 2 keV Chandra (blue)
contours overlaid for the H$\beta$, [O III]$\lambda\lambda 4960,5008$, [O
I]$\lambda 6300$, H$\alpha$, [N II]$\lambda\lambda 6548,6584$, and [S
II]$\lambda\lambda 6718,6733$ emission lines. MUSE images have a pixel size of
0.4 arcsec/pixel. 8.4 GHz VLA contours were drawn at 5, 20, and 50 times the
rms level of the background. Chandra contours were smoothed with a 1.5″
Gaussian kernel radius and drawn in blue at 0.25$\,\cdot\,$10-15 erg
cm${}^{-2}\,$s-1. Flux maps have units of 10-17 erg cm${}^{-2}\,$s-1, while
velocity and velocity dispersion maps have units of km s-1.


# Charge Density Wave Surface Reconstruction
in van der Waals Layered Materials

Sung-Hoon Lee<EMAIL_ADDRESS>Department of Applied Physics and Institute of
Natural Sciences, Kyung Hee University, Yongin 17104, Korea  Doohee Cho
<EMAIL_ADDRESS>Department of Physics, Yonsei University, Seoul 03722,
Republic of Korea

###### Abstract

Surface reconstruction plays a vital role in determining the surface
electronic structure and chemistry of semiconductors and metal oxides.
However, it has been commonly believed that surface reconstruction does not
occur in van der Waals layered materials, as they do not undergo significant
bond breaking during surface formation. In this study, we present evidence
that charge density wave (CDW) order in these materials can, in fact, cause
CDW surface reconstruction through interlayer coupling. Using density
functional theory calculations on the 1T-TaS2 surface, we reveal that CDW
reconstruction, involving collective and modest atomic displacements in the
subsurface layer, results in a significant modification of the surface
electronic structure, transforming it from a Mott insulator to a band
insulator. This new form of surface reconstruction explains several previously
unexplained observations on the 1T-TaS2 surface and has important implications
for interpreting surface phenomena in CDW-ordered layered materials.

##  1 Introduction

The transition metal dichalcogenide 1T-TaS2 has been an exclusive material
platform for investigating the charge density wave (CDW) order and strong
electron correlation in two dimensions (2D) 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11.
Recently, It has become the focus of attention due to the competition between
in-plane electron correlation and interlayer coupling of CDWs 12, 13, 14, 15,
16, 17, 18, 19, 20, 21, 22, 23, 24. The triangular lattice of Ta ions in the
$d^{1}$ electronic configuration that makes up a 1T-TaS2 layer, triggers CDW
instabilities. When cooled to below $\sim$180 K, the Star of David (SoD) CDW
distortion with a $(\sqrt{13}\times\sqrt{13})$ periodicity is observed. While
the twelve Ta ions in the enlarged unit cell form a CDW gap, the remaining Ta
ion at the center of the SoD creates a half-filled state. These half-filled
localized states in a single layer result in a 2D Mott insulator 2, 25, 26
with a peak-to-peak gap of $\sim$0.5 eV between the Hubbard bands 27, 28. Bulk
1T-TaS2 is also an insulator with a bandgap of $\sim$0.1 eV 29 and was
previously thought to be a Mott insulator with negligible interlayer coupling
2, 3, 4, 5, 6, 7, 8, 9, 10, 11. However, recent density functional theory
(DFT) calculations 12, 13, 14, 15 show that interlayer coupling is
substantial, and CDW stacking along the out-of-plane direction has a
significant impact on the electronic structure. One particular bilayer CDW
stacking configuration, called AL stacking, has been identified as the most
stable 15 and has been confirmed by various experimental studies 30, 31, 32,
33, 16, 17, 18. This AL stacking causes interlayer dimerization between the
localized states, making the system a band insulator 14, 15.

Studies of the 1T-TaS2 surface have uncovered a more intricate interaction
between electron correlation and CDW stacking. Scanning tunneling microscopy
(STM) studies 19, 20, 21 have identified three different surfaces: two
insulating surfaces with peak-to-peak gaps of 0.4 eV and 0.24 eV, and metallic
surfaces. The difference was attributed to variations in CDW stacking on the
surface 19. By considering that bulk CDW stacking can have two cleavage
planes, the bilayer-terminated surface was attributed to the large-gap
insulating surface, and the single-layer-terminated surface was attributed to
the small-gap insulating surface. Although a many-body theoretical study
provided support for these assignments 22, a subsequent STM study 23 raised an
issue with this understanding. Both the upper and lower terraces of a single-
layer surface step exhibited the same large gap, despite one of the two
terraces having to be single-layer-terminated in the case of bilayer CDW
stacking. To reconcile this inconsistency, it was suggested that the large gap
is also a single-layer feature and a correlation-driven Mott gap 23. On the
other hand, angle-resolved photoemission spectroscopy (ARPES) data 8, 34
showed two bands near the Fermi level with distinct $k_{z}$ dispersions: a
dispersive band due to interlayer coupling 35, 15 and a flat band at around
$-0.2$ eV. Given that the peak at $-0.2$ eV originates from a large-gap
surface featuring peaks at $\pm 0.2$ eV 19, 20, 21, the flatness of this peak
along the $k_{z}$ direction remained unexplained.

Here, using density functional calculations, we investigate the interplay
between electron correlation and interlayer coupling on the surface of CDW-
bilayer-stacked 1T-TaS2. Our investigation reveals a previously unknown
surface reconstruction, which explains the puzzling observations on the
surface. Specifically, we unequivocally determine the CDW structures of the
large-gap, small-gap, and metallic surfaces, including their respective energy
stabilities and spin structures. The large-gap surface is characterized by a
band-insulating bilayer-terminated surface, while the small-gap surface is a
Mott-insulating single-layer-terminated surface, which interestingly is less
stable than the former. Consequently, when the area of the less stable small-
gap surface is large, a spontaneous CDW reconstruction takes place, with the
Mott-insulating single layer shifting into the third layer and leaving a CDW
bilayer on the surface. The resulting two types of bilayer-terminated
surfaces, one with a subsurface 2D Mott insulator and the other without,
explain the inconsistent STM observation and flat ARPES band. In addition, we
identify metallic surfaces as Mott insulators that can be described by
trilayer Hubbard models. Overall, our results provide a solid foundation for
understanding the complex interaction between electron correlation and
interlayer coupling in 1T-TaS2. The CDW surface reconstruction we propose
involves moderate atomic displacements, yet it causes significant changes in
the surface electronic structure and is highly likely to occur in other CDW-
ordered layered materials as well.

##  2 Results

In this study, we adopt the CDW stacking notation from Ref. 15. Due to the
three-fold rotational symmetry of the SoD cluster (Fig. 1a), it has five
distinct stacking interfaces denoted by stacking vectors
$\mathbf{T}_{s}=\mathbf{c}$, $\mathbf{a}+\mathbf{c}$,
$2\mathbf{a}+\mathbf{c}$, $-2\mathbf{a}+\mathbf{c}$, and
$-\mathbf{a}+\mathbf{c}$. These interfaces are labeled as A, B, C, L, and M
(refer to Fig. 1c,d) and are written in italics. It is worth noting that the C
and L interfaces are different when considering the S sublattice, which is not
displayed in Fig. 1c.

Fig. 1:  Star of David (SoD) distortion and CDW stacking order in 1T-TaS2. a
Atomic structure of 1T-TaS2 under SoD distortion. b Isodensity surface of the
half-filled state that is localized at the center of the SoD cluster. c
Labeling convention for Ta atoms in an SoD cluster. d Arrangement of half-
filled states for A, L, and AL stacking, with markings indicating the location
of the SoD’s center for each layer on the left, and the CDW stacking interface
between the layers on the right.  Fig. 2:  Stability of bulk 1T-TaS2 as a
function of CDW stacking and the Hubbard $U$ potential. The black lines
represent the results obtained from nonmagnetic (NM) calculations, while the
blue lines represent the results obtained from calculations with layer-by-
layer antiferromagnetic (AFM) order.

###  2.1 Determination of $U$ for surface calculations

Before investigating surface properties, we evaluate the stability of bulk
1T-TaS2 in relation to CDW stacking and the Hubbard $U$ potential to determine
the appropriate range of $U$ values in our DFT + $U$ calculations (Fig. 2).
For nonmagnetic cases, L stacking is found to be the most stable among the
five single-layer stacking configurations for all considered $U$ values.
Furthermore, bilayer AL stacking is even more stable than L stacking for all
$U$ values. This confirms and extends previous findings (with $U=0$ eV) that
the ground-state CDW configuration for 1T-TaS2 is AL stacking when nonmagnetic
15. The high stability of the L interface originates from two factors:
enhanced van der Waals energy gain and favorable interlayer S-S bonding, with
AL stacking further increasing stability through a bandgap opening and
accompanying electronic energy gain 15. When magnetic order is included, L
stacking begins to favor layer-by-layer antiferromagnetic (AFM) order at
$U\approx 1.4$ eV and becomes the most stable stacking configuration above
$U_{c}=1.92$ eV. The onset of layer-by-layer AFM order for A and AL stacking
is delayed to $U\approx 2.5$ and $2.1$ eV respectively, due to stronger
interlayer coupling at the A interface. Therefore, the ground-state CDW
stacking phase is nonmagnetic AL stacking for $U<U_{c}$ and layer-by-layer
AFM-ordered L stacking for $U>U_{c}$. The experimental evidence for bilayer AL
stacking 30, 31, 32, 16, 17, 18, 33 sets a limit of $U<U_{c}$ in our DFT + $U$
calculations. Our calculations suggest that the $U$ values of 2.27–2.94 eV
calculated in previous DFT + $U$ studies 12, 36, 37 (obtained using Dudarev et
al.’s scheme 38 for Ta $5d$ orbitals in 1T-TaS2) are too high and may
overestimate electron correlation. We use a moderate $U$ value of 1.25 eV in
our subsequent surface calculations.

Fig. 3:  Electronic structures of the large-gap and small-gap insulating
surfaces. a Nonmagnetic band structure of an eight-layer slab with AL
stacking. b Projected density of states (PDOS) of the central Ta $5d_{z^{2}}$
orbital in the top layer (left) and the STM $dI/dV$ spectrum 19 (right). c
Nonmagnetic band structure (left), spin-polarized band structure (center), and
spin density (right) of a seven-layer slab with AL stacking. The majority and
minority spin bands are represented by red and green lines, respectively, in
the spin-polarized band structure. The shaded area (cyan, orange, and violet)
associated with each band indicates the spectral weight of the central Ta
$5d_{z^{2}}$ orbital in the top layer.

###  2.2 Large-gap and small-gap insulating surfaces

We begin by investigating two basic surfaces of 1T-TaS2: one cleaved between
bilayers and the other within a bilayer. We model these surfaces using eight-
and seven-layer slabs, respectively, both with a ground-state AL stacking. The
eight-layer slab, with bilayer termination on both surfaces, is a band
insulator that lacks any magnetic order; it does not exhibit magnetic order,
even with $U=1.75$ eV. The bandgap of the slab is determined by the layers
inside, and it is calculated to be 0.06 eV, which is nearly equivalent to the
bulk bandgap of 0.05 eV. The spectral weight of the Ta $5d_{z^{2}}$ orbital at
the SoD’s center in the top layer, as shown by the shaded area in cyan in Fig.
3a, indicates that the surface bilayer has a slightly increased interlayer-
hybridization gap due to a minor surface relaxation. The corresponding
projected density of states (PDOS) manifests two peaks at $-0.18$ and $0.08$
eV, exhibiting a peak-to-peak gap of 0.26 eV. Given that DFT calculations
typically underestimate the bandgap substantially, the peak positions and
intensities of the calculated PDOS match well with the STM $dI/dV$ spectrum of
a large-gap (0.4 eV) insulating surface (Fig. 3b, blue lines). This finding
implies that the peaks at $\pm 0.2$ eV identified by STM and ARPES correspond
to the bonding and antibonding hybridization bands of the surface bilayer.

The seven-layer slab with an undimerized top layer exhibits a narrow half-
filled band in the nonmagnetic calculation (Fig. 3c, left panel). When spin
polarization is included, the half-filled band becomes fully gapped, causing
the split Hubbard bands to merge with the bulk bands (Fig. 3c, center panel).
The calculated PDOS of the Ta $5d_{z^{2}}$ orbital in the top layer displays a
peak-to-peak gap of 0.16 eV between the upper and lower Hubbard bands and a
double-peak structure in the occupied states. These PDOS features are well-
represented in the STM spectrum of a small-gap (0.24 eV) insulating surface
(Fig. 3b, red lines). The spatial spin distribution (Fig. 3c, right panel)
indicates that 88% of the spin is confined to the top layer. These results
indicate that the single CDW layer on the surface is a 2D Mott insulator that
is decoupled from the bulk.

###  2.3 Correlated surface trilayers and “metallic” surfaces

Fig. 4:  Electronic structures of the “metallic” surfaces. a, c Spin-polarized
band structure (left), PDOS (center), and spin density (right) of the single-
layer surface with the A interface a and the M interface c. b, d Nonmagnetic
band structure (left) and PDOS (right) of the single-layer surface with the A
interface b and the M interface d. In the spin-polarized band structures, the
majority and minority spin bands are represented by red and green lines,
respectively. The shaded area (orange, violet, and cyan) associated with each
band indicates the spectral weight of the central Ta $5d_{z^{2}}$ orbital in
the top layer.

The extent of decoupling between the surface single layer and subsurface
layers depends on the interlayer hopping constant ($t_{\perp}$), which is
determined by the stacking interface between them. The value of $t_{\perp}$ is
the largest at the A interface ($\mathbf{T}_{s}=\mathbf{c}$), while it is the
smallest at the L or C interface ($\mathbf{T}_{s}=\mp
2\mathbf{a}+\mathbf{c}$). It is instructive to consider the strongly coupled
case of A-interface termination, wherein a surface trilayer is formed (Fig.
4a). Even in this case, the band structure suggests that the surface acts as a
Mott insulator. Notably, the spin of the occupied Hubbard band in this case is
divided equally between the first and third layers, bypassing the second
layer. This behavior can be explained by regarding the surface trilayer as a
2D lattice of vertical trimers with small lateral hopping constants. A single
trimer with equal nearest-neighbor hopping integrals of $t$ has three
eigenstates with eigenvalues of $\pm\sqrt{2}t$ and 0, and the zero-energy
eigenstate is $(1,0,-1)$ (see Supplementary Note 1). The half-filled zero-
energy eigenstates of the vertical trimers form the half-filled band of the
surface trilayer, as evidenced by the nonmagnetic band structure and PDOS in
Fig. 4b. Furthermore, the nonmagnetic band structure illustrates that the
energy level corresponding to $-\sqrt{2}t$ is $-0.22$ eV, indicating that the
interlayer hopping constant at the A interface, $t_{\perp}^{A}$, is $0.16$ eV.
It also shows that the bandwidth of the half-filled band is only 0.07 eV,
implying that the in-plane hopping constant ($t_{\parallel}$) is much smaller
than $t_{\perp}^{A}$. When spin polarization is switched on, the half-filled
band splits into upper and lower Hubbard bands with a peak-to-peak gap of 0.07
eV (Fig. 4a). The large decrease in the gap compared to that of the single-
layer surface with the L interface (0.16 eV) is caused by the doubled spatial
extent of the half-filled states.

The experimental studies using STM have suggested that the single-layer
surfaces with the M or B interface ($\mathbf{T}_{s}=\mp\mathbf{a}+\mathbf{c}$)
exhibit “metallic” spectra 19, 23. The surface electronic structure of the
M-interface case (Fig. 4c) lies between those of the L\- and A-interface
cases: the three surface layers form a 2D correlated system, with
asymmetrically split spin in the first (71%) and third (26%) layers of the
surface trilayer. The trimer model with asymmetric nearest-neighbor hopping
integrals $t_{\perp}^{M}$ and $t_{\perp}^{A}$ yields an estimate of 0.61 for
the ratio $t_{\perp}^{M}/t_{\perp}^{A}$ (Supplementary Note 1), indicating
that $t_{\perp}^{M}$ is around 0.10 eV. The calculated top-layer PDOS shows
good agreement with the STM $dI/dV$ spectrum of a “metallic” surface, as shown
in Fig. S1. It is important to note that the split Hubbard bands of the
surface trilayer, with a peak-to-peak gap of 0.09 eV, and the absence of a
metallic peak at the Fermi level in the STM spectrum indicate that the
“metallic” surface is not a metal but instead a Mott insulator. The single-
layer surfaces with the B and C interfaces exhibit similar characteristics to
those with the M and L interfaces, respectively (Fig. S2).

###  2.4 Surface energetics

In Fig. 5a, we compare the surface energies of the different surfaces we have
analyzed so far. The surfaces are divided into three groups based on their
stability: the bilayer-terminated surface is the most stable, while the
single-layer surface of the L-interface is unstable, and all other single-
layer surfaces, including the trilayer (A-interface) case, are even more
unstable. These results are consistent with the energetics for CDW stacking in
the bulk, as shown in Fig. 2, and provide a rationale for observations from
STM studies. Previous research 7, 8, 9, 10, 19, 20, 21, 23 has found that most
surfaces, regardless of whether they were cooled to or cleaved at cryogenic
temperatures, exhibit a large gap, which is consistent with the bilayer-
terminated surface being the most stable. The small-gap insulating surface,
identified as the L-interface single-layer surface, has only been observed in
small domains near step edges or domain boundaries 19, 23, unless stabilized
by alkali metal adatoms 21, due to its instability. Finally, “metallic”
surfaces are less frequent than the small-gap surface 23, reflecting their
high instability. Overall, the calculated surface energies provide valuable
insights into the stability and features of the various surfaces.

Fig. 5:  Surface energies and a buried single-layer surface. a Surface
energies of various surfaces [meV$/(\sqrt{13}\times\sqrt{13})$]. b Spin-
polarized band structure (left), PDOS (center), and spin density (right) of a
buried single-layer surface. The majority and minority spin bands are
represented by red and green lines, respectively, in the spin-polarized band
structure. The shaded area in cyan (orange and violet) on the band structure
represents the spectral weight of the central Ta $5d_{z^{2}}$ orbital in the
first (third) layer.

###  2.5 CDW surface reconstruction

The STM observations of a large gap in both the upper and lower terraces of
single-layer steps 23 can be explained by the presence of a buried single
layer. If a single CDW layer on the surface is covered by a CDW bilayer with L
stacking interface (Fig. 5b), then the surface energy is lowered by 32
meV$/(\sqrt{13}\times\sqrt{13})$, making it comparable to the energy of a
plain bilayer-terminated surface without a subsurface single layer (Fig. 5a).
Although the buried single layer forms a 2D Mott insulator in the third layer
with a peak-to-peak gap of 0.12 eV, the PDOS on the surface bilayer is almost
indistinguishable from that of the plain bilayer-terminated surface, with a
peak-to-peak gap of 0.25 eV. This indicates that STM measurements would
observe the same large gap of 0.4 eV, with peaks at $\pm 0.2$ eV, on bilayer-
terminated surfaces, regardless of the presence of a single CDW layer in the
third layer. The STM observations can be explained as follows: when a single-
layer step is created by cleaving an AL-stacked bulk, one of the two terraces
becomes single-layer-terminated (Fig. 6a, upper panel). If the size of the
single-layer-terminated terrace is large enough, the CDWs undergo spontaneous
reconstruction, causing the single CDW layer in the first layer to move to the
third layer (Fig. 6a, lower panel). The final stabilized step has bilayer
termination on both terraces (refer to Fig. S3, for the CDW reconstruction in
two consecutive steps).

Fig. 6:  CDW surface reconstruction. a Surface reconstruction of the
L-interface single-layer surface on the right side of a single-layer step into
a buried single-layer surface. It results in a reduction in the surface energy
of 32 meV$/(\sqrt{13}\times\sqrt{13})$. b CDW shift in the second layer. It
induces a displacement of the center of the SoD from the ‘A’ site (blue SoD)
to the ‘L’ site (red SoD), resulting in a shift of 5.8 Å. The magnified image
of the region encircled in green illustrates the lateral displacement of a Ta
atom by 0.4 Å. The gray triangular lattice represents the undistorted lattice
points of the Ta sublattice.

The proposed CDW reconstruction on the surface of 1T-TaS2 exhibits two
remarkable features. Firstly, it involves a CDW shift only in the second
layer, while the other layers remain unaltered (Fig. 6a). Despite this, the
surface electronic structure undergoes a drastic transformation from a Mott
insulator to a band insulator. Secondly, the CDW shift occurs through a
lateral translation of the SoD rather than a uniform atomic displacement (Fig.
6b). This process involves a collection of modest displacements of atoms
around their undistorted lattice points, with a maximum lateral displacement
of 0.4Å, yet results in a 5.8 Å shift in the SoD position. This collective
atomic motion without addition or subtraction of atoms is a hallmark of CDW
surface reconstruction and contrasts with the typical surface reconstructions
observed on semiconductor surfaces, which involve adding or subtracting atoms
to satisfy local charge states and minimize the number of broken bonds 39.

##  3 Discussion

The identification of CDW reconstruction in single-layer steps represents an
important finding. It reveals that the majority of the 1T-TaS2 surface, which
exhibits a large gap, is not simply the plain bilayer-terminated surface, but
rather comprises two types of bilayer-terminated surfaces: the plain one and
the one with a subsurface single CDW layer. This finding holds significant
implications for surface observations, not only with STM but also with all
surface-sensitive tools.

The observation of two bands with distinct $k_{z}$ dispersions around $-0.2$
eV in the ARPES measurements 8, 34 can be attributed directly to the presence
of two types of bilayer-terminated surfaces. Our findings indicate that the
plain bilayer-terminated surface is the origin of the dispersive band, while
the flat band arises from the bilayer surface with a subsurface single layer.
Because of Coulomb repulsion, the Mott-insulating subsurface layer prevents
electronic hopping between the surface bilayer and the subsurface bulk,
causing the surface bilayer to be decoupled from the bulk and exhibit a flat
dispersion along the $k_{z}$ direction. It is important to note that the flat
band at $-0.2$ eV cannot be attributed to the Mott-insulating surface single
layer because the corresponding STM $dI/dV$ spectrum shows a broad double-peak
structure, as depicted in Fig. 3b.

Our findings also shed new light on surface observations made through time-
resolved photoemission spectroscopy (tr-PES) and scanning thermoelectric
microscopy (SThEM). The previous tr-PES measurements 11 observed an
ultrashort-lived peak around $+0.2$ eV and broad spectral features in the gap
region. A recent theoretical study utilizing dynamical mean-field theory
simulations 40 attributed the sharp peak to the single-layer-terminated
surface and the in-gap states to the plain bilayer-terminated surface,
assuming that the former constituted a substantial portion of the surface.
However, our results demonstrate that this assumption is incorrect, and we
suggest that the peak at $+0.2$ eV originates from a bilayer surface that is
decoupled from the bulk. The presence of peaks at $\pm 0.2$ eV in the bilayer
surface corroborates this interpretation. Similarly, in the SThEM study 24,
unexpected variations in thermopower were ascribed to the Mott-insulating
single-layer-terminated surface, but our findings suggest that this
interpretation needs to be reevaluated.

Lastly, we note that the 1T-TaS2 surface presents a promising opportunity to
study the trilayer Hubbard model. To understand the impact of interlayer
coupling on the strongly correlated 2D system of cuprate high-temperature
superconductors, the bilayer Hubbard model has been extensively studied, and
the model at half-filling has been found to undergo a transition from a Mott
insulator to a band insulator as interlayer coupling increases 41, 42, 43.
However, the trilayer Hubbard model, which more accurately represents the
properties of multiply stacked materials in three dimensions, has received
less attention. Our investigation of the 1T-TaS2 surface has revealed that the
“metallic” surfaces, which are in fact Mott-insulating, correspond to trilayer
Hubbard model systems in the strong interlayer coupling regime ($t_{\perp}\gg
t_{\parallel}$). These findings indicate that the 1T-TaS2 surface is a rare
system for exploring the properties of the trilayer Hubbard model, and inspire
further research in this area.

In conclusion, our investigation of the surface electronic properties of
1T-TaS2 has led to the discovery of a previously unknown phenomenon of surface
reconstruction of CDWs in van der Waals layered materials, with significant
implications for the interpretation of experimental observations using various
surface probes. This finding suggests that CDW surface reconstruction may be
widespread in other CDW-ordered layered materials, making it an area worthy of
further investigation. Additionally, our theoretical work has revealed a
fascinating and intricate landscape of strongly correlated phenomena in
1T-TaS2, highlighting its potential as an excellent model system for studying
the physics of correlated layered materials, especially the interplay between
CDW, electron correlation, and interlayer coupling.

##  4 Methods

###  4.1 Density functional theory calculations

We performed DFT + $U$ calculations using the Vienna ab initio simulation
package (VASP) 44, 45. These calculations employed the projector-augmented-
wave method 46, the generalized gradient approximation 47, and the DFT + $U$
scheme developed by Dudarev et al 38. To account for van der Waals
interactions, we used the Tkatchenko-Scheffler approach 48. The electronic
wave functions were expanded using a plane-wave basis set, with cutoff
energies of 323 eV and 259 eV for bulk and surface calculations, respectively.
We carried out k-space integration using a $4\times 4\times 8$ mesh in the
Brillouin zone of the $13\times 13\times 1$ supercell. To optimize the lattice
constants and atomic positions for each stacking configuration, we used the
variable cell optimization method in VASP. For AL stacking and $U=1.25$ eV for
Ta $5d$ orbitals, the calculated lattice constants were $a=3.36$ Å and
$c=5.77$ Å, compared to experimental values of $a=3.37$ Å and $c=5.90$ Å18.
The agreement between theory and experiment is quite good, especially given
that the $c$ parameter is mostly determined by van der Waals interactions. To
model the surface, we employed a periodic slab geometry, where each slab
comprised seven or eight TaS2 layers, and the thickness of the vacuum region
exceeded $14$ Å. We relaxed all atoms until all residual forces were less than
0.01 eV/Å.

## References

  * Wilson et al. 1975 Wilson, J.; Di Salvo, F.; Mahajan, S. Charge-density waves and superlattices in the metallic layered transition metal dichalcogenides. _Adv. Phys._ 1975, _24_ , 117–201
  * Fazekas and Tosatti 1979 Fazekas, P.; Tosatti, E. Electrical, structural and magnetic properties of pure and doped 1T-TaS2. _Philos. Mag. Part B_ 1979, _39_ , 229–244
  * Rossnagel 2011 Rossnagel, K. On the origin of charge-density waves in select layered transition-metal dichalcogenides. _J. Phys.: Condens. Matter_ 2011, _23_ , 213001
  * Sipos et al. 2008 Sipos, B.; Kusmartseva, a. F.; Akrap, a.; Berger, H.; Forró, L.; Tutis, E. From Mott state to superconductivity in 1T-TaS2. _Nat. Mater._ 2008, _7_ , 960–965
  * 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–180
  * Yu et al. 2015 Yu, Y.; Yang, F.; Lu, X. F.; Yan, Y. J.; Cho, Y.-H.; Ma, L.; Niu, X.; Kim, S.; Son, Y.-W.; Feng, D.; Li, S.; Cheong, S.-W.; Chen, X. H.; Zhang, Y. Gate-tunable phase transitions in thin flakes of 1T-TaS2. _Nat. Nanotechnol._ 2015, _10_ , 270–276
  * Kim et al. 1994 Kim, J.-J.; Yamaguchi, W.; Hasegawa, T.; Kitazawa, K. Observation of Mott localization gap using low temperature scanning tunneling spectroscopy in commensurate 1T-TaS2. _Phys. Rev. Lett._ 1994, _73_ , 2103–2106
  * Ma et al. 2016 Ma, L.; Ye, C.; Yu, Y.; Lu, X. F.; Niu, X.; Kim, S.; Feng, D.; Tománek, D.; Son, Y.-W.; Chen, X. H.; Zhang, Y. A metallic mosaic phase and the origin of Mott-insulating state in 1T-TaS2. _Nat. Commun._ 2016, _7_ , 10956
  * Cho et al. 2016 Cho, D.; Cheon, S.; Kim, K.-S.; Lee, S.-H.; Cho, Y.-H.; Cheong, S.-W.; Yeom, H. W. Nanoscale manipulation of the Mott insulating state coupled to charge order in 1T-TaS2. _Nat. Commun._ 2016, _7_ , 10453
  * Qiao et al. 2017 Qiao, S.; Li, X.; Wang, N.; Ruan, W.; Ye, C.; Cai, P.; Hao, Z.; Yao, H.; Chen, X.; Wu, J.; Wang, Y.; Liu, Z. Mottness collapse in 1T-TaS2-xSex tansition-metal dichalcogenide: An interplay between localized and itinerant orbitals. _Phys. Rev. X_ 2017, _7_ , 041054
  * Ligges et al. 2018 Ligges, M.; Avigo, I.; Golež, D.; Strand, H.; Beyazit, Y.; Hanff, K.; Diekmann, F.; Stojchevska, L.; Kalläne, M.; Zhou, P.; Rossnagel, K.; Eckstein, M.; Werner, P.; Bovensiepen, U. Ultrafast doublon dynamics in photoexcited 1T-TaS2. _Phys. Rev. Lett._ 2018, _120_ , 166401
  * Darancet et al. 2014 Darancet, P.; Millis, A. J.; Marianetti, C. A. Three-dimensional metallic and two-dimensional insulating behavior in octahedral tantalum dichalcogenides. _Phys. Rev. B_ 2014, _90_ , 045134
  * Ritschel et al. 2015 Ritschel, T.; Trinckauf, J.; Koepernik, K.; Büchner, B.; Zimmermann, M. v.; Berger, H.; Joe, Y. I.; Abbamonte, P.; Geck, J. Orbital textures and charge density waves in transition metal dichalcogenides. _Nat. Phys._ 2015, _11_ , 328–331
  * Ritschel et al. 2018 Ritschel, T.; Berger, H.; Geck, J. Stacking-driven gap formation in layered 1T-TaS2. _Phys. Rev. B_ 2018, _98_ , 195134
  * Lee et al. 2019 Lee, S.-H.; Goh, J. S.; Cho, D. Origin of the insulating phase and first-order metal-insulator transition in 1T-TaS2. _Phys. Rev. Lett._ 2019, _122_ , 106404
  * von Witte et al. 2019 von Witte, G.; Kißlinger, T.; Horstmann, J. G.; Rossnagel, K.; Schneider, M. A.; Ropers, C.; Hammer, L. Surface structure and stacking of the commensurate $(\sqrt{13}\times{}\sqrt{13})$R$13.9^{\circ}$ charge density wave phase of 1T-TaS${}_{2}(0001)$. _Phys. Rev. B_ 2019, _100_ , 155407
  * Stahl et al. 2020 Stahl, Q.; Kusch, M.; Heinsch, F.; Garbarino, G.; Kretzschmar, N.; Hanff, K.; Rossnagel, K.; Geck, J.; Ritschel, T. Collapse of layer dimerization in the photo-induced hidden state of 1T-TaS2. _Nat. Commun._ 2020, _11_ , 1247
  * Wang et al. 2020 Wang, Y. D.; Yao, W. L.; Xin, Z. M.; Han, T. T.; Wang, Z. G.; Chen, L.; Cai, C.; Li, Y.; Zhang, Y. Band insulator to Mott insulator transition in 1T-TaS2. _Nat. Commun._ 2020, _11_ , 4215
  * Butler et al. 2020 Butler, C. J.; Yoshida, M.; Hanaguri, T.; Iwasa, Y. Mottness versus unit-cell doubling as the driver of the insulating state in 1T-TaS2. _Nat. Commun._ 2020, _11_ , 2477
  * Butler et al. 2021 Butler, C. J.; Yoshida, M.; Hanaguri, T.; Iwasa, Y. Doublonlike excitations and their phononic coupling in a Mott charge-density-wave system. _Phys. Rev. X_ 2021, _11_ , 011059
  * Lee et al. 2021 Lee, J.; Jin, K. H.; Yeom, H. W. Distinguishing a Mott insulator from a trivial insulator with atomic adsorbates. _Phys. Rev. Lett._ 2021, _126_ , 196405
  * Petocchi et al. 2022 Petocchi, F.; Nicholson, C. W.; Salzmann, B.; Pasquier, D.; Yazyev, O. V.; Monney, C.; Werner, P. Mott versus hybridization gap in the low-temperature phase of 1T-TaS2. _Phys. Rev. Lett._ 2022, _129_ , 016402
  * Wu et al. 2022 Wu, Z.; Bu, K.; Zhang, W.; Fei, Y.; Zheng, Y.; Gao, J.; Luo, X.; Liu, Z.; Sun, Y.-P.; Yin, Y. Effect of stacking order on the electronic state of 1T-TaS2. _Phys. Rev. B_ 2022, _105_ , 035109
  * Kim et al. 2022 Kim, D.; Shin, E.-C.; Lee, Y.; Lee, Y. H.; Zhao, M.; Kim, Y.-H.; Yang, H. Atomic-scale thermopower in charge density wave states. _Nat. Commun._ 2022, _13_ , 4516
  * Yi et al. 2018 Yi, S.; Zhang, Z.; Cho, J.-H. Coupling of charge, lattice, orbital, and spin degrees of freedom in charge density waves in 1T-TaS2. _Phys. Rev. B_ 2018, _97_ , 041413(R)
  * Pasquier and Yazyev 2022 Pasquier, D.; Yazyev, O. V. Ab initio theory of magnetism in two-dimensional 1T-TaS2. _Phys. Rev. B_ 2022, _105_ , L081106
  * Lin et al. 2020 Lin, H.; Huang, W.; Zhao, K.; Qiao, S.; Liu, Z.; Wu, J.; Chen, X.; Ji, S.-H. Scanning tunneling spectroscopic study of monolayer 1T-TaS2 and 1T-TaSe2. _Nano Res._ 2020, _13_ , 133–137
  * Vaňo et al. 2021 Vaňo, V.; Amini, M.; Ganguli, S. C.; Chen, G.; Lado, J. L.; Kezilebieke, S.; Liljeroth, P. Artificial heavy fermions in a van der Waals heterostructure. _Nature_ 2021, _599_ , 582–586
  * Gasparov et al. 2002 Gasparov, L. V.; Brown, K. G.; Wint, A. C.; Tanner, D. B.; Berger, H.; Margaritondo, G.; Gaál, R.; Forró, L. Phonon anomaly at the charge ordering transition in 1T-TaS2. _Phys. Rev. B_ 2002, _66_ , 094301
  * Nakanishi and Shiba 1984 Nakanishi, K.; Shiba, H. Theory of three-dimensional orderings of charge-density waves in 1T-TaX2 (X: S, Se). _J. Phys. Soc. Jpn._ 1984, _53_ , 1103–1113
  * Naito et al. 1986 Naito, M.; Nishihara, H.; Tanaka, S. Nuclear magnetic resonance and nuclear quadrupole resonance study of 181Ta in the commensurate charge density wave state of 1T-TaS2. _J. Phys. Soc. Jpn._ 1986, _55_ , 2410–2421
  * Ganal et al. 1990 Ganal, P.; Butz, T.; Lerf, A.; Naito, M.; Nishihara, H. The 181Ta nuclear quadrupole interaction in the charge density wave phases of 1T-TaS2. _Z. Naturforsch. A_ 1990, _45_ , 439–444
  * Ishiguro and Sato 1991 Ishiguro, T.; Sato, H. Electron microscopy of phase transformations in 1T-TaS2. _Phys. Rev. B_ 1991, _44_ , 2046–2060
  * Jung et al. 2022 Jung, J.; Park, J. W.; Kim, J.; Yeom, H. W. Surface enhanced electron correlation on the trivial quasi-two-dimensional bulk insulator 1T-TaS2. _Phys. Rev. B_ 2022, _106_ , 155406
  * Ngankeu et al. 2017 Ngankeu, A. S.; Mahatha, S. K.; Guilloy, K.; Bianchi, M.; Sanders, C. E.; Hanff, K.; Rossnagel, K.; Miwa, J. A.; Breth Nielsen, C.; Bremholm, M.; Hofmann, P. Quasi-one-dimensional metallic band dispersion in the commensurate charge density wave of 1T-TaS2. _Phys. Rev. B_ 2017, _96_ , 1–5
  * Martino et al. 2020 Martino, E.; Pisoni, A.; Ćirić, L.; Arakcheeva, A.; Berger, H.; Akrap, A.; Putzke, C.; Moll, P. J.; Batistić, I.; Tutiš, E.; Forró, L.; Semeniuk, K. Preferential out-of-plane conduction and quasi-one-dimensional electronic states in layered 1T-TaS2. _npj 2D Mater. Appl._ 2020, _4_ , 1–9
  * Boix-Constant et al. 2021 Boix-Constant, C.; Mañas-Valero, S.; Córdoba, R.; Baldoví, J. J.; Rubio, A.; Coronado, E. Out-of-plane transport of 1T-TaS2/graphene-based van der Waals heterostructures. _ACS Nano_ 2021, _15_ , 11898–11907
  * Dudarev et al. 1998 Dudarev, S. L.; Botton, G. A.; Savrasov, S. Y.; Humphreys, C. J.; Sutton, A. P. Electron-energy-loss spectra and the structural stability of nickel oxide: An LSDA+U study. _Phys. Rev. B_ 1998, _57_ , 1505–1509
  * Zangwill 1988 Zangwill, A. _Physics at surfaces_ ; Cambridge University Press, 1988
  * Petocchi et al. 2022 Petocchi, F.; Chen, J.; Li, J.; Eckstein, M.; Werner, P. Photo-induced charge dynamics in 1T-TaS2. _arXiv:2212.09583_ 2022,
  * Gall et al. 2021 Gall, M.; Wurz, N.; Samland, J.; Chan, C. F.; Köhl, M. Competing magnetic orders in a bilayer Hubbard model with ultracold atoms. _Nature_ 2021, _589_ , 40–43
  * Fuhrmann et al. 2006 Fuhrmann, A.; Heilmann, D.; Monien, H. From Mott insulator to band insulator: A dynamical mean-field theory study. _Phys. Rev. B_ 2006, _73_ , 245118
  * Kancharla and Okamoto 2007 Kancharla, S. S.; Okamoto, S. Band insulator to Mott insulator transition in a bilayer Hubbard model. _Phys. Rev. B_ 2007, _75_ , 193103
  * Kresse and Furthmüller 1996 Kresse, G.; Furthmüller, J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. _Phys. Rev. B_ 1996, _54_ , 11169–11186
  * Kresse and Joubert 1999 Kresse, G.; Joubert, D. From ultrasoft pseudopotentials to the projector augmented-wave method. _Phys. Rev. B_ 1999, _59_ , 1758–1775
  * Blöchl 1994 Blöchl, P. E. Projector augmented-wave method. _Phys. Rev. B_ 1994, _50_ , 17953–17979
  * Perdew et al. 1996 Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized gradient approximation made simple. _Phys. Rev. Lett._ 1996, _77_ , 3865–3868
  * Tkatchenko and Scheffler 2009 Tkatchenko, A.; Scheffler, M. Accurate molecular van der Waals interactions from ground-state electron density and free-atom reference data. _Phys. Rev. Lett._ 2009, _102_ , 073005

### Supplementary Note 1. Tight-binding model of a trimer

In a tight-binding framework, the Hamiltonian of a trimer with nearest-
neighbor hopping integrals, $t_{s}$ and $t_{b}$, can be expressed as

$H=\begin{pmatrix}0&t_{s}&0\\\ t_{s}&0&t_{b}\\\ 0&t_{b}&0\end{pmatrix},$ (1)

and its three eigenvectors and eigenvalues are

$\left|E_{\pm}=\pm\sqrt{t_{s}^{2}+t_{b}^{2}}\right\rangle=\begin{pmatrix}t_{s}\\\
\mp\sqrt{t_{s}^{2}+t_{b}^{2}}\\\
t_{b}\end{pmatrix},\qquad\left|E_{0}=0\right\rangle=\begin{pmatrix}t_{b}\\\
0\\\ -t_{s}\end{pmatrix}.$ (2)

The solutions for three special cases are

  * •

$t_{s}/t_{b}\approx 1$ (corresponding to the A-interface single-layer
(trilayer) surface)

$\left|E_{\pm}=\pm\sqrt{2}t_{b}\right\rangle=\begin{pmatrix}1\\\
\mp\sqrt{2}\\\
1\end{pmatrix},\qquad\left|E_{0}=0\right\rangle=\begin{pmatrix}1\\\ 0\\\
-1\end{pmatrix}.$ (3)

  * •

$t_{s}/t_{b}\approx 3/5$ (corresponding to the M-interface single-layer
surface)

$\left|E_{\pm}=\pm\frac{\sqrt{34}}{5}t_{b}=\pm
1.17t_{b}\right\rangle=\begin{pmatrix}3/5\\\ \mp\sqrt{34}/5\\\
1\end{pmatrix},\qquad\left|E_{0}=0\right\rangle=\begin{pmatrix}1\\\ 0\\\
-3/5\end{pmatrix}.$ (4)

  * •

$t_{s}/t_{b}\approx 1/3$ (corresponding to the L-interface single-layer
surface)

$\left|E_{\pm}=\pm\frac{\sqrt{10}}{3}t_{b}=\pm
1.05t_{b}\right\rangle=\begin{pmatrix}1/3\\\ \mp\sqrt{10}/3\\\
1\end{pmatrix},\qquad\left|E_{0}=0\right\rangle=\begin{pmatrix}1\\\ 0\\\
-1/3\end{pmatrix}.$ (5)

Fig. S1:  Comparison of a the STM $dI/dV$ data (from Nat. Commun. 11, 2477
(2020)) and b the corresponding PDOS data of this work. It was
undistinguishable whether the stacking interface of Region 2 was M or B. Here,
we compare Region 2 with the M-interface single-layer surface.

Fig. S2:  Band structure, PDOS, and spin density of the single-layer surface
for four stacking interfaces: a B, b M, c C, and d L. In the band structure,
the majority and minority spin bands are depicted by red and green lines,
respectively. The shaded area (orange and violet) associated with each band
represents the spectral weight of the central Ta $5d_{z^{2}}$ orbital in the
top layer.  Fig. S3:  CDW surface reconstruction for two distinct types of
single-layer steps (top panels) and for two consecutive steps consisting of
both types (bottom panel). The CDW reconstruction leads to a reduction in the
surface energy by 32 meV$/(\sqrt{13}\times\sqrt{13})$. This process entails a
CDW shift in the second layer without any modifications to the other layers.