# SOFIA FEEDBACK Survey: The Pillars of Creation in [C II] and Molecular Lines

Ramsey L. Karim University of Maryland, Department of Astronomy, Room 1113 PSC
Bldg. 415, College Park, MD 20742-2421, USA  Marc W. Pound University of
Maryland, Department of Astronomy, Room 1113 PSC Bldg. 415, College Park, MD
20742-2421, USA  Alexander G.G.M. Tielens University of Maryland, Department
of Astronomy, Room 1113 PSC Bldg. 415, College Park, MD 20742-2421, USA
Maitraiyee Tiwari University of Maryland, Department of Astronomy, Room 1113
PSC Bldg. 415, College Park, MD 20742-2421, USA  Max-Planck Institute for
Radioastronomy, Auf dem Hügel, D-53121 Bonn, Germany  Lars Bonne SOFIA Science
Center, USRA, NASA Ames Research Center, M.S. N232-12, Moffett Field, CA
94035, USA  Mark G. Wolfire University of Maryland, Department of Astronomy,
Room 1113 PSC Bldg. 415, College Park, MD 20742-2421, USA  Nicola Schneider I.
Physikalisches Institut, Universität zu Köln, Zülpicher Str. 77, D-50937 Köln,
Germany  Ümit Kavak SOFIA Science Center, USRA, NASA Ames Research Center,
M.S. N232-12, Moffett Field, CA 94035, USA  Lee G. Mundy University of
Maryland, Department of Astronomy, Room 1113 PSC Bldg. 415, College Park, MD
20742-2421, USA  Robert Simon I. Physikalisches Institut, Universität zu Köln,
Zülpicher Str. 77, D-50937 Köln, Germany  Rolf Güsten Max-Planck Institute for
Radioastronomy, Auf dem Hügel, D-53121 Bonn, Germany  Jürgen Stutzki I.
Physikalisches Institut, Universität zu Köln, Zülpicher Str. 77, D-50937 Köln,
Germany  Friedrich Wyrowski Max-Planck Institute for Radioastronomy, Auf dem
Hügel, D-53121 Bonn, Germany  Netty Honingh  I. Physikalisches Institut,
Universität zu Köln, Zülpicher Str. 77, D-50937 Köln, Germany

###### Abstract

We investigate the physical structure and conditions of photodissociation
regions (PDRs) and molecular gas within the Pillars of Creation in the Eagle
Nebula using SOFIA FEEDBACK observations of the [C II] 158 $\mu$m line. These
observations are velocity resolved to 0.5 ${\rm km~{}s^{-1}}$ and are analyzed
alongside a collection of complimentary data with similar spatial and spectral
resolution: the [O I] 63 $\mu$m line, also observed with SOFIA, and rotational
lines of CO, HCN, HCO+, CS, and N2H+. Using the superb spectral resolution of
SOFIA, APEX, CARMA, and BIMA, we reveal the relationships between the warm PDR
and cool molecular gas layers in context of the Pillars’ kinematic structure.

We assemble a geometric picture of the Pillars and their surroundings informed
by illumination patterns and kinematic relationships and derive physical
conditions in the PDRs associated with the Pillars. We estimate an average
molecular gas density $n_{{{\rm H_{2}}}}\sim{1.3\times 10^{5}}$ cm-3 and an
average atomic gas density $n_{{\rm H}}\sim{1.8\times 10^{4}}$ cm-3 and infer
that the ionized, atomic, and molecular phases are in pressure equilibrium if
the atomic gas is magnetically supported. We find pillar masses of 103, 78,
103, and 18 $M_{\odot}$ for P1a, P1b, P2, and P3 respectively, and evaporation
times of $\sim$1–2 Myr. The dense clumps at the tops of the pillars are
currently supported by the magnetic field. Our analysis suggests that
ambipolar diffusion is rapid and these clumps are likely to collapse within
their photoevaporation timescales.

††facilities: SOFIA(upGREAT), APEX(LAsMA, SEPIA660), CARMA, BIMA,
JWST(NIRCam), Herschel(PACS, SPIRE), Spitzer(IRAC)††software: astropy (Astropy
Collaboration et al., 2013, 2018, 2022), Spectral Cube
(https://doi.org/10.5281/zenodo.3558614), regions
(https://doi.org/10.5281/zenodo.7259631), pvextractor
(https://github.com/radio-astro-tools/pvextractor), numpy (Harris et al.,
2020), scipy (Virtanen et al., 2020), pandas (McKinney, 2010, 2011),
matplotlib (Hunter, 2007), pdrtpy (Pound & Wolfire, 2023)

##  1 Introduction

The Eagle Nebula, also known as M16, is a well-studied H II region lying above
the Galactic midplane at $(l,~{}b)=(16\arcdeg 9540,~{}+0\arcdeg 7934)$,
$(\alpha,\,\delta)_{J2000}=(18\mbox{${}^{\mathrm{h}}$}18\mbox{${}^{\mathrm{m}}$}48\fs
0,~{}$-$13\arcdeg 48\arcmin 24\arcsec)$. The H II region is illuminated by a
$\sim$2 Myr old stellar cluster, NGC 6611 (Hillenbrand et al., 1993), born
from the giant molecular cloud (GMC) W 37 (Zhan et al., 2016). Filaments
primarily extended in the Galactic $b$ direction, parallel to the larger scale
elongation of W 37, lie within W 37 near the location of M16 (Hill et al.,
2012; Xu et al., 2019).

The H II region associated with M 16 is powered by the early-type stars in the
NGC 6611 cluster, a $\sim 2\times 10^{4}~{}M_{\odot}$ cluster whose most
massive member is estimated to be $\sim 80~{}M_{\odot}$ (Hillenbrand et al.,
1993; Hill et al., 2012). The cluster, at a heliocentric distance of $\sim$2
kpc, comprises a number of pre-main sequence (PMS) stars as well as stars
surrounded by circumstellar material, both signs of recent or ongoing star
formation. NGC 6611 members exhibit some age spread, as evidenced by several
hundred PMS stars whose ages range from 0.25–1 Myr and a massive main sequence
population about 2 Myr old lying alongside a handful of evolved massive stars,
at least one of which is B2.5 I, which is around 6 Myr old (Hillenbrand et
al., 1993). Guarcello et al. (2007, 2009, 2010) and De Marchi et al. (2013)
delve into the star formation history of the region and characterize two
different populations of stars: a group of $>$8 Myr old stars which likely
existed before the formation of the M16 H II region, and a group of $\sim$1
Myr-old stars, broadly consistent with the $\sim$2 Myr age determination of
NGC 6611 by other authors, whose formation event may have been triggered by a
supernova shell’s arrival $\sim$3 Myr ago. We adopt the distance of $1740\pm
130$ pc determined using GAIA parallax observations (Kuhn et al., 2019) and a
main sequence age of 2 Myr based on studies from the last two decades (Belikov
et al., 1999, 2000; Dufton et al., 2006; Martayan et al., 2008; Stoop et al.,
2023).

Within M16 lies an iconic pillar system whose Hubble Space Telescope (HST),
and now JWST, images are well known to the public (Hester et al., 1996). The
three main Pillars, extending towards a handful of bright O5–7 stars $\sim 3$
pc away in projection (Hillenbrand et al., 1993), are seen roughly parallel to
the plane of the sky as in Figure 1. Spectroscopic studies by Pound (1998) and
McLeod et al. (2015) have concluded that all the Pillars are inclined slightly
towards or away from the observer, and the tallest pillar, called Pillar 1
(P1), is actually a superposition of two pillars: the northern half, called
P1a, is actually behind the cluster and the southern half, called P1b, is in
between the cluster and the observer (see Figure 1). Pound (1998) detects
coherent molecular gas flows along the line of sight, which they conclude are
projected motions along the pillar. They find that the observed radial
velocity gradient along the pillar implies that the dynamical timescale for
dissipation, $\sim 10^{5}$ yr, may be shorter than their estimated
photoevaporation timescale $\sim 10^{7}$ yr. Both Pound (1998) and McLeod et
al. (2015) suggest that Pillars 2 and 3 (P2 and P3) are between the cluster
and the observer and that P2 points away from the observer towards the cluster
(throughout this paper, “points” refers to the orientation of the head; i.e.,
the head of P2 is farther from the observer than the tail). The line-of-sight
(LOS) velocity data collected by Pound (1998) and McLeod et al. (2015)
indicate that P3 points towards the observer, and the illumination data from
the HST observations presented by Hester et al. (1996) and the MUSE
observations analyzed by McLeod et al. (2015) are consistent with P3 lying
between the cluster and the observer. Together, the velocity and illumination
data imply that P3 does not point towards the cluster. P3 might point towards
a cluster member displaced towards the observer along the line of sight
(Section 4.2).

Due to the high spatial resolution of the JWST images, structures within the
Pillar system are well-defined. We use the JWST images throughout this study
to measure the on-sky angular widths of features.

Figure 1: (Left) JWST color composite prepared by Joseph DePasquale (STScI),
Anton M. Koekemoer (STScI), Alyssa Pagan (STScI). Colors are Purple: F090W,
Blue: F187N, Cyan: F200W, Yellow: F335M, Orange: F444W, Red: F470N. The image
was obtained from
https://webbtelescope.org/contents/media/images/2022/052/01GF423GBQSK6ANC89NTFJW8VM.
(Center) Three-color composite using JWST filters F090W (blue), F187N (green),
and F335M (red). The stretches are nonlinear and limits have been adjusted.
The 3.3 µm band, in red, includes a PAH feature and therefore indicates far-
ultraviolet illumination. Section 2.4 includes a more detailed summary of the
significance of these filters. (Right) Schematic diagram of the Pillars on the
sky at the same angular scale as the two images to the left. Features are
marked with labels which will be used throughout the paper. P1a refers to the
northern half of Pillar 1, including the Cap and Eastern and Western Threads.
P1b refers to the southern half and includes the Eastern and Western Horns and
part of the Shared Base. The dashed line to the south marks a boundary which
is kinematically discontinuous with the Ridge despite its apparent continuity
in the visible and IR images; this is discussed in Section 3.6. The three
images are nearly RA-Dec aligned; “up” in the image is 2.3∘ east of north. The
scale bar in the top-left corner shows $1\arcmin\approx 0.5$ pc at 1740 pc
line of sight.

###  1.1 The Photodissociation Region

The photodissociation region (PDR) is the region of far-ultraviolet (FUV;
6–13.6 eV) illuminated neutral atomic and molecular gas just behind the
ionization front (see Tielens & Hollenbach 1985 and Wolfire et al. 2022 for
detailed background). Due primarily to photoelectric heating from small grains
and polycyclic aromatic hydrocarbons (PAHs) (Bakes & Tielens, 1994), the
atomic region at the cloud surface is warmer ($T\sim 100$ K) than the
molecular gas ($T\sim 10$ K) found deeper in the cloud. Most of the cooling in
the atomic PDR is done through a small handful of fine-structure lines,
principally the 158 µm [C II] line (singly-ionized carbon), and at $n\gtrsim
10^{4}$ cm-3, the 63 µm and 145 µm [O I] lines (neutral oxygen) (Tielens &
Hollenbach, 1985; Hollenbach & Tielens, 1999; Wolfire et al., 2022).
Spectroscopically resolved observations of these lines, such as those we
present here, allow us to probe the kinematics of the warmer, outer atomic
region of the PDR. Observations of [C II] and [O I] are presented by Schneider
et al. (2012, 2021) towards pillars and globules in Cygnus X.

At $A_{V}\gtrsim 5$, carbon is found primarily in the molecular phase as
carbon monoxide (CO). From $A_{V}\sim 5$–10, FUV radiation still warms the
predominantly molecular gas, which we can probe with low- to mid-J CO lines.
These lines, and particularly the bright 12CO(J=1-0) line, should be optically
thick and trace this molecular PDR layer as opposed to the colder, less
illuminated molecular gas $A_{V}\gtrsim 10$.

Molecular lines more sensitive to dense gas ($n\gtrsim 10^{5}$ cm-3) will
probe the denser, colder layers beyond the PDR which are not heated by FUV
radiation. By comparing observations of dense gas tracers, such as the HCO+,
HCN, and CS observed towards the Pillars by White et al. (1999), with the PDR
tracers introduced above, we can explore not only the molecular gas inventory
of the Pillars but also how mass moves between phases and leaves the Pillars
through bulk flows (Pound, 1998) or photoevaporation (Hester et al., 1996;
McLeod et al., 2015).

We present a multi-wavelength analysis using both velocity-resolved and
continuum observations tracing a variety of gas phases within and around the
Pillars, from the cold, dense layers deep within the pillar heads to the warm,
outer layers illuminated by the bright members of NGC 6611. This study
presents the first velocity-resolved [C II] and [O I] line observations of the
Pillars of Creation which probe the conditions and kinematics of the FUV-
illuminated PDR layer between the ionization front and the molecular gas
within the Pillars. We describe the Pillars and surrounding features in all
tracers in Section 3 and determine the location of major PDRs and discuss the
illumination geometry in Section 4. Our derivation of column densities, number
densities, and pressures in the atomic and molecular phases of the gas is
discussed in Section 5, where we also discuss pressure equilibrium between
these phases and the ionized gas. We discuss the photoevaporative timescale of
the Pillars in Section 6 and include a summary of our work and some closing
remarks in Section 7.

##  2 Observations

###  2.1 SOFIA

M16 was observed with upGREAT111upGREAT and GREAT were developed by the MPI
für Radioastronomie and the KOSMA/Universität zu Köln, in cooperation with the
MPI für Sonnensystemforschung and the DLR Institut für Planetenforschung. on
SOFIA on 9 flights between 2019 and 2022 from Palmdale, California and Tahiti
in the [C II] 3P${}_{3/2}\rightarrow$ 3P1/2 transition at 158 µm and in the [O
I] 63 µm line in parallel with the upGREAT receiver (Risacher et al., 2018)
onboard SOFIA. Note that we did not use the large undersampled [O I] map from
the FEEDBACK mapping but instead, data from an earlier PI program (see below).
An area of $\sim 590$ arcmin2 was mapped in the on-the-fly (OTF) mode and
atmospheric calibration was done with the GREAT pipeline (Guan et al., 2012).
A Fast Fourier Transform Spectrometer (FFTS) with 4 GHz instantaneous
bandwidth and a frequency resolution of 0.244 MHz serves as a backend (Klein
et al., 2012).

The nominal angular resolution of the [C II] and [O I] data is 14.1″ and 6″,
respectively, but here we use a [C II] data cube with a spatial resolution of
15.4″, a grid of 3.5″, and a spectral binning of 0.5 ${\rm km~{}s^{-1}}$. The
noise RMS in one channel is typically 1.0 K (Table 1). A 3rd order baseline
was removed from the spectra, which were then averaged with 1/2 weighting
(baseline noise). Spectra are presented on a main beam brightness temperature
scale $T_{\rm MB}$ using an average main beam efficiency of 0.65. The forward
efficiency is $\eta_{\rm f}$ = 0.97. See Schneider et al. (2020) for more
observational details. These [C II] observations are made public through the
NASA/IPAC Infrared Science Archive
(IRSA)222https://irsa.ipac.caltech.edu/Missions/sofia.html.

A smaller area of M16, covering mostly the Pillars of Creation, was observed
during 2 flights in October 2016 (Cycle 5) with the GREAT receiver with two
channels. The 7-pixel HFA was tuned to the [O I] 63 $\mu$m line and the single
pixel L2 channel to the CO(J=16-15) line at 1841.345506 GHz, which we do not
discuss in this study. We employed the Fast Fourier Transform Spectrometer
backend AFFTS. The center IF (intermediate frequency) was 1455 MHz for the [O
I] line and 1000 MHz for the CO line. The map was obtained in beam-switched
on-the-fly mapping mode. The stepsize of the map was 2.4″ which is at a higher
sampling than the FEEDBACK large mapping. The angular resolution of the [O I]
data is 6″. All line intensities are reported as main beam temperatures scaled
with main-beam efficiencies of 0.69 and 0.68 for [O I] and CO, respectively,
and a forward efficiency of 0.97. From the spectra, a 3rd order baseline was
removed and spectra were then averaged with 1/2 weighting (baseline noise). We
smoothed the [O I] data to a resolution of 0.4 ${\rm km~{}s^{-1}}$.

###  2.2 CARMA and BIMA

We used the Combined Array for Research in Millimeter-wave Astronomy (CARMA,
Bock et al., 2006) to map the Eagle Nebula pillars in four spectral lines that
trace high H2 volume density gas. The observations were obtained with the
15-element array comprised of the six 10.4 m antennas and nine 6.1 m antennas.
Three E-array configuration tracks were obtained on August 8, September 8, and
September 10, 2012, and one D-array configuration track was obtained on
November 12, 2012. The antenna signals were transmitted to the eight-band
spectral line correlator. Four bands were used to observe the spectral lines
HCN(J=1-0), HCO+(J=1-0), N2H+(J=1-0), and CS(J=2-1) in 7.8 MHz bandwidths with
spectral resolution 24 kHz/channel ($\Delta$V $\sim$0.08 ${\rm km~{}s^{-1}}$).
Four bands were used to measure continuum at $\sim 92$ GHz, each with 490 MHz
bandwidth and 12.5 MHz/channel. Combining the two sidebands gave a total
continuum bandwidth of 3.84 GHz. (Continuum emission was detected at the level
of $\sim$ 10 mJy; it is not discussed in this study.) The Eagle was observed
with a 37-point hexagonal mosaic centered on
$(\alpha,\,\delta)_{J2000}=(18\mbox{${}^{\mathrm{h}}$}18\mbox{${}^{\mathrm{m}}$}51\fs
29,-13\arcdeg 15\arcmin 02.32\arcsec)$. The map was sampled at the Nyquist
interval of the 10.4 m antennas at the CS rest frequency because it is the
highest frequency and thus smallest interval. Finally, we also use in this
paper the CO(J=1-0), 13CO(J=1-0) and C18O(J=1-0) archival Berkeley-Illinois-
Maryland Array (BIMA) data of the Eagle pillars from Pound (1998) and Pound et
al. (2007).

The data were reduced using the MIRIAD package (Sault et al., 1995, 2011).
After phase, amplitude, passband, and flux calibration, and flagging of bad
data, visibilities were inverted onto a 0.5″ spatial grid and 0.1 ${\rm
km~{}s^{-1}}$ channels using a robust weighting value of zero (Briggs, 1995).
The inverted images were deconvolved with the MIRIAD task mossdi which uses
CLEAN algorithm of Steer et al. (1984). Deconvolved maps were restored with a
fitted 2D Gaussian beam. Details of the observations are provided in Table 1.

The N2H+(J=1-0) line splits into a series of hyperfine transitions, the two
strongest of which, (J, F1, F) $=$ (1–0, 2–1, 2–1) and (1–0, 2–1, 3–2), lie
$\sim$1 ${\rm km~{}s^{-1}}$ from each other. Since these are not well
separated given the $\sim$1 ${\rm km~{}s^{-1}}$ full-width at half-max (FWHM)
of the lines, we only use spatial information from them and do not use their
velocities. Our observations include a transition of N2H+ (J, F1, F) $=$ (1–0,
0–1, 1–2) which is well-separated from other transitions given the typical
FWHM, so we use velocity information from this line. The HCN(J=1-0) line is
also split into hyperfine transitions, but the strong central transition is
well separated from the others given the typical FWHM.

###  2.3 APEX

M16 was mapped on September 18-20, 2019, in good weather conditions
(precipitable water vapor pwv = 0.6–0.9 mm) in the 13CO(J=3-2) and 12CO(J=3-2)
transitions using the LAsMA spectrometer on the APEX333APEX, the Atacama
Pathfinder Experiment is a collaboration between the Max-Planck-Institut für
Radioastronomie, Onsala Space Observatory (OSO), and the European Southern
Observatory (ESO). telescope (Güsten et al., 2006). LAsMA is a 7-pixel single
polarization heterodyne array that allows simultaneous observations of the two
isotopomers in the upper (12CO) and lower (13CO) sideband of the receiver,
respectively.

The array is arranged in a hexagonal configuration around a central pixel with
a spacing of about two beam widths ($\theta_{\rm MB}=18.2\arcsec$ at 345.8
GHz) between the pixels. It uses a K mirror as de-rotator. The backends are
advanced Fast Fourier Transform Spectrometers (Klein et al., 2012) with a
bandwidth of $2\times 4$ GHz and a native spectral resolution of 61 kHz. The
mapping was done in total power on-the-fly mode using a reference position at
$(\alpha,\,\delta)_{J2000}=(18\mbox{${}^{\mathrm{h}}$}20\mbox{${}^{\mathrm{m}}$}46\fs
3,-13\arcdeg 14\arcmin 56\arcsec)$.

The mapped region, centered at
$18\mbox{${}^{\mathrm{h}}$}18\mbox{${}^{\mathrm{m}}$}35\fs 7,-13\arcdeg
43\arcmin 31.0\arcsec$ and of size $30\arcmin\times 22\arcmin$ at $-40$ deg
angle (CCW against positive RA), was split into $2\times 2$ tiles. Each tile
was scanned with a spacing of 9″ (oversampling to 6″ in scanning direction),
resulting in a uniformly sampled map with high fidelity. All spectra are
calibrated in TMB (main-beam efficiency $\eta_{\rm MB}=0.68$ at 345.8 GHz).
All linear baselines were removed and all data re-sampled to 0.1 ${\rm
km~{}s^{-1}}$ spectral bins. The final data cubes are built with a pixel size
of 9.1″ and an 18.2″ beam after gridding.

M16 was mapped on May 3 and June 28, 2021, in reasonable weather conditions
(pwv = 0.5–0.9 mm) in the 12CO(J=6-5) transition using the SEPIA660 receiver
and the FFTS1 backend spectrometer on the APEX telescope. The instrument was
tuned so that the 12CO(J=6-5) line rest frequency 691.473 GHz lay in the upper
sideband, whose bandwidth is 4 GHz and frequency resolution is 61 kHz. The
OTF-mapped region is a $4.3\arcmin\times 4.3\arcmin$ square centered at
$18\mbox{${}^{\mathrm{h}}$}18\mbox{${}^{\mathrm{m}}$}51\fs 2,-13\arcdeg
50\arcmin 0\arcsec$. Spectra are calibrated in TMB (main-beam efficiency
$\eta_{\rm MB}=0.45$). The final cube is gridded to 4.5″ pixels with a 9.6″
beam and has a channel width of 0.25 ${\rm km~{}s^{-1}}$.

###  2.4 Ancillary Data

M16 was observed at 8 µm using the InfraRed Array Camera (IRAC, Fazio et al.
2004) on board Spitzer as part of the GLIMPSE program (Benjamin et al., 2003),
and at 70 µm and 160 µm using the Photodetector Array Camera and Spectrometer
(PACS, Poglitsch et al. 2010) and 250, 350, and 500 µm using the Spectral and
Photometric Imaging REceiver (SPIRE, Griffin et al. 2010) aboard the Herschel
Space Observatory (Pilbratt et al., 2010) as part of the Hi-GAL Galactic plane
survey (Molinari et al., 2010). The IRAC image was obtained from the GLIMPSE
website444https://irsa.ipac.caltech.edu/data/SPITZER/GLIMPSE/ and the PACS and
SPIRE observations were obtained from the Herschel Science Archive (HSA), and
they are accessible through the NASA/IPAC IRSA: GLIMPSE Team (2020) and Hi-GAL
Team (2020a, b, c, d, e). Ultraviolet (UV) radiation, emitted in abundance
from NGC 6611, excites large hydrocarbon molecules called polycyclic aromatic
hydrocarbons, or PAHs, which fluoresce in the infrared (IR), and the 8 µm
filter covers a particularly strong feature in the PAH spectrum (Tielens,
2008). Detailed studies by Indebetouw et al. (2007) and Flagey et al. (2011)
discuss these and other Spitzer mid-IR images. Far-infrared (FIR) 70 and 160
µm emission traces warm dust illuminated with FUV radiation from the stars,
with 70 µm relatively more sensitive to temperature than column density; M16
was studied in detail in the FIR by Hill et al. (2012).

The Pillars were observed using JWST NIRCam (Rieke et al., 2023) as part of
the Cycle 1 outreach campaign (PI: Pontoppidan, PID #2739) and their data made
publicly available via the Mikulski Archive for Space Telescopes (MAST):
10.17909/fbc0-1930 (catalog 10.17909/fbc0-1930). We use images in the filters
F090W, F187N, and F335M, which trace 0.9 µm continuum and background emission
from starlight, the 1.87 µm Pa-$\alpha$ recombination line from ionized
hydrogen, and the 3.3 µm PAH feature. The 0.9 and 1.87 µm observations are
particularly useful for locating areas of high near-IR (NIR) extinction
towards the Pillars, and the 3.3 µm observation highlights illuminated PDR
surfaces at $\sim 200$ AU resolution.

Table 1: Observations summary.

Species | Frequency | Beam Size | Beam PA | dV | RMS
---|---|---|---|---|---
| (GHz) | (″) | (degrees) | (${\rm km~{}s^{-1}}$) | (K)
Continuum | 92.0000000 | $9.25\times 5.94$ | –5.34 |  $\cdots$ | 0.0042
HCN | 88.6318470 | $11.31\times 7.22$ | –6.57 | 0.10 | 0.52
HCO+ | 89.1885180 | $11.30\times 7.15$ | –7.30 | 0.10 | 0.56
N2H+ | 93.1735050 | $10.59\times 6.76$ | –6.48 | 0.10 | 0.64
CS | 97.9809680 | $10.12\times 6.45$ | –6.29 | 0.10 | 0.73
C18O(J=1-0) | 109.7821600 | $13.73\times 7.42$ | –9.55 | 0.267 | 0.60
13CO(J=1-0) | 110.2013530 | $6.80\times 4.29$ | –8.15 | 0.266 | 2.06
12CO(J=1-0) | 115.2712040 | $7.51\times 4.39$ | –1.38 | 0.254 | 5.03
13CO(J=3-2) | 330.5879653 | 20.0 | 0 | 0.111 | 0.68
12CO(J=3-2) | 345.7959899 | 19.2 | 0 | 0.106 | 0.55
12CO(J=6-5) | 691.4730000 | 9.6 | 0 | 0.25 | 1.65
12CO(J=16-15) | 1841.345506 | 15.9 | 0 |  $\cdots$ |  $\cdots$
[C II] | 1900.536900 | 15.4 | 0 | 0.5 | 1.0
[O I] | 4744.777490 | 6.7 | 0 | 0.4 | 1.9

Note. — Beam PA is the position angle of the elliptical beam, measured in
degrees east of north. dV is the velocity bin width used in this study. RMS is
the root-mean-squared noise within each velocity bin for the given beam and
dV. The 92 GHz continuum and CO(J=16-15) line were observed but are not used
in this study.

##  3 Results

We present a rich collection of observations of the Pillars of Creation
alongside publicly available archival data. Integrated intensity and continuum
maps in Figure 2 show the Pillars in a variety of tracers, from warm atomic
gas ([C II] and [O I]) and UV illuminated gas (8 µm PAH and 70 µm dust
continuum) to warm (CO) and dense (HCO+ and CS) molecular gas.

This data set is ripe for comparison between tracers of different physical
conditions as well as between our velocity-resolved line maps and extremely
high spatial resolution JWST images. Comparing integrated [C II] intensity to
the F335M image in Figure 3 reveals myriad relationships between emission
features both strong and weak, from the bright emission along the Pillars to
the faint legs extending from the base of the Pillars down to the southeast
(near
$\alpha,\,\delta=18\mbox{${}^{\mathrm{h}}$}18\mbox{${}^{\mathrm{m}}$}56\fs,\,-13\arcdeg
51\arcmin 20\arcsec$ in Figure 3), while the CS integrated intensity map
overlaid in that figure unveils the density structure beneath.

Pillar 1 appears between ${\rm V_{LSR}}=24$–26 ${\rm km~{}s^{-1}}$ in the
channel maps in Figure 4 and Pillars 2 and 3 appear between ${\rm
V_{LSR}}=20$–23 ${\rm km~{}s^{-1}}$. The $\sim$10 pc scale filaments observed
by Hill et al. (2012) and Xu et al. (2019) lie at similar velocities ${\rm
V_{LSR}}=20$, 22.5, and 25 ${\rm km~{}s^{-1}}$. Much of the [C II] emission in
the channel maps is diffuse emission surrounding the Pillars between ${\rm
V_{LSR}}\sim 24$–27 ${\rm km~{}s^{-1}}$, and so the diffuse component line
overlaps with and contaminates the P1 component line. We describe in Appendix
B the case and procedure for subtracting a background from [C II] spectra
towards the Pillars to handle this contamination. The diffuse 3.3 µm haze,
apparent by comparing the upper left corner of Figure 3 to the slightly dimmer
lower right corner, may trace the same diffuse PDR as this ${\rm V_{LSR}}\sim
25$ ${\rm km~{}s^{-1}}$ [C II] component. The pillar system connects towards
the southeast to a $\sim$10 pc scale system of emission features (only a small
part of which is shown in Figure 4) at ${\rm V_{LSR}}\sim 25~{}{\,{\rm
km~{}s^{-1}}}$ which may be the illuminated edge of the H II region. At the
lower velocities of P2 and P3, a different set of emission features extend
towards the southeast, away from the stars, highlighted in blue in Figure 5.
Higher velocity channel maps ${\rm V_{LSR}}\geq 29~{}{\,{\rm km~{}s^{-1}}}$
reveal a faint, ring-like feature $\sim$2 pc in diameter with a bright
southern edge. Only the western side of this ring is included in the Figure 4
channel maps. Since it has a significantly different morphology than the
pillar system, we deem it to be a foreground or background feature unrelated
to the Pillars but likely part of the H II region and a candidate for future
work.

The [C II] and [O I] line profiles (FWHM $\sim 2$–3 ${\rm km~{}s^{-1}}$) match
the expected line width for atomic gas PDR emission. The peak main beam
temperatures $T_{\rm MB}$ of [C II] reach $\sim$40 K towards the Pillars. The
12CO line widths $\sim 2$ ${\rm km~{}s^{-1}}$ and peak $T_{\rm MB}$ (reaching
$\sim$80 K in the (J=1-0) line) indicate their origin in the warm molecular
gas in the PDR. The HCO+, HCN, and CS line profiles are relatively thinner
(FWHM $\sim 1$–2 ${\rm km~{}s^{-1}}$) and reach lower $T_{\rm MB}$ ($\sim$20 K
in HCO+ and HCN; $\sim$10 K in CS). All observed line widths are likely
dominated by turbulent motions. No evidence of self-absorption in [C II] or [O
I] is detected towards the bright Shared Base or any other location towards
the Pillars (see Guevara et al. 2020 for a discussion of [C II] self-
absorption in other Galactic star-forming regions).

Throughout the paper, we refer to Pillars 1, 2, and 3 as P1, P2 and P3. As was
discussed in Section 1, P1 refers to the entire structure of Pillar 1
including both P1a and P1b (McLeod et al., 2015). The “head” of a pillar
refers to the region of the pillar nearer to the illuminating stars, and the
“tail” or “base” refers to the region farther from the stars. The “body”
refers to the bulk midsection of the pillar. “Along” the body or “parallel” to
the Pillars refers to the radial direction with respect to the stars, while
“across” the body or “transverse” refers to the tangential direction with
respect to the stars. We use the proper noun “Pillars” to refer to the system
in M16 which is the subject of this work and the common noun “pillars” to
refer to the generic structure found in other regions as well.

Figure 2: Integrated intensity and photometry centered on the Pillars.
Contoured observations are marked with (c). The beam for each observation is
shown in the lower right corner; the beam for the contoured observation, if
present, is shown to the left of the color-scale observation’s beam. All
observations are shown at the same angular scale. [C II] and [O I] are
integrated between ${\rm V_{LSR}}=18$–27 ${\rm km~{}s^{-1}}$, as they tend to
have longer blue wings, and molecular lines not otherwise specified are
integrated between ${\rm V_{LSR}}=19$–27 ${\rm km~{}s^{-1}}$. N2H+ and HCN are
integrated between ${\rm V_{LSR}}=12.6$–32 ${\rm km~{}s^{-1}}$ with respect to
the rest frequencies given in Table 1 in order to include several satellite
lines for each species (see Section 2). Contours are [0.5, 3.0] K ${\rm
km~{}s^{-1}}$ for C18O(J=1-0); [4, 24, 44] K ${\rm km~{}s^{-1}}$ for
13CO(J=3-2); and [9, 18, 27, 45] $\times 10^{5}$ MJy sr-1 for the 160 µm image
(the last contour spacing is intentionally uneven to increase visibility
towards P1a). The contour colors are chosen to increase contrast with the
image and have no further significance. The dashed white contour on the
interferometric molecular line observations marks 50% gain from the primary
beam mosaic pattern; outside this contour the sensitivity falls off (noise
increases). Note the BIMA CO(J=1-0) maps are smaller than the CARMA maps of
the other species and so cover the southeast Ridge with very low sensitivity.
Nonetheless it is weakly detected. Figure 3: Integrated [C II] (black) and CS
(white) line intensities between ${\rm V_{LSR}}=19$–27 ${\rm km~{}s^{-1}}$
overlaid on the F335M image. [C II] integrated intensities are marked at the
[25, 50, 75, 100, 125, 150, 175, 200, 225] K ${\rm km~{}s^{-1}}$ levels, and
CS at [0.8, 3.2, 5.6, 8, 10.4, 12.8, 15.2, 17.6, 20, 22.4, 24.8] K ${\rm
km~{}s^{-1}}$. The dotted white line marks the 50% gain contour for the CS
observations; see the caption of Figure 2 for more detail. The [C II] beam is
shown in grey in the lower right corner, and the CS beam superimposed in
white. Note the spatial offsets between the two lines towards the Cap and the
Threads, among other locations. The large-scale [C II] emission indicated by
the lowest two contours is the ${\rm V_{LSR}}\sim 25~{}{\,{\rm km~{}s^{-1}}}$
background discussed in detail in Appendix B. Figure 4: 158 µm [C II] line
channel maps binned to 1 ${\rm km~{}s^{-1}}$ and centered on the Pillars. The
beam is shown in the lower right corner of each map. Figure 5: (Left) Color
composite using [C II] line integrated intensities between ${\rm
V_{LSR}}=19$–21.5 ${\rm km~{}s^{-1}}$ (blue), 22–23.5 ${\rm km~{}s^{-1}}$
(green), and 24–27.5 ${\rm km~{}s^{-1}}$ (red). The color stretches are linear
and start at 0 but have different upper limits. (Right) NOAO optical composite
prepared by T.A. Rector (NRAO/AUI/NSF and NOIRLab/NSF/AURA) and B.A. Wolpa
(NOIRLab/NSF/AURA) using observations from the WIYN 0.9 m telescope at the
Kitt Peak National Observatory. The image was obtained from
https://noirlab.edu/public/images/noao-04086 and coordinate metadata applied
using Astrometry.neta (Lang et al., 2010). The colors show the O III line at
499 nm (blue), the H$\alpha$ line at 656 nm (green), and the S II line at 672
nm (red). ahttps://astrometry.net Figure 6: Spectra of the [C II], [O I],
CO(J=6-5), and CS lines towards the three brightest [O I] positions. The
lower-right three panels show the absolute spectra in Kelvins. The upper-left
three panels show the same spectra normalized to their maxima. Half-maximum,
for estimating the FWHM, is marked with dashed horizontal lines. The [C II]
spectra are shown after the spectral background subtraction described in
Appendix B, and the unsubtracted [C II] spectra are shown with dotted lines in
the lower-right panels for context. Vertical lines mark every 1 ${\rm
km~{}s^{-1}}$ between 20–27 ${\rm km~{}s^{-1}}$. The [O I], CO, and CS
observations are all convolved to the [C II] beam, shown in the lower right
corner of the central panel. The central panel shows [C II] (color) and [O I]
(contour) emission integrated between 20–27 ${\rm km~{}s^{-1}}$. The [O I]
contours are placed at [0, 6, 12, 18, 24, 30] K ${\rm km~{}s^{-1}}$. Dashed
grey lines in the central panel show the area covered by these [O I]
observations. Velocities are relative to the local standard of rest (LSR).

###  3.1 P1

The [C II] and [O I] lines, all observed molecular lines, and the 3–500 µm
continuum peak in brightness twice along the length of Pillar 1: towards P1a
and P1b. The peak brightness ratio of P1a to P1b is larger for denser gas
tracers and smallest for [C II], in which P1b is brighter. We show spectra of
[C II], [O I], CO(J=6-5), and CS towards these brightness peaks in Figure 6.
Between the peaks, warm tracers like [C II] and 3–160 µm continuum are
continuous and remain brighter than the background along the pillar, whereas
the molecular lines trace the discontinuous, clumpy structure observed by
Thompson et al. (2002) in near-IR images and White et al. (1999) in their
radio, sub-mm, and IR observations (compare the [C II] and CS in Figure 3).

A bright rim of PAH and optical emission lies atop the molecular emission
towards P1a due to direct illumination from NGC 6611 over a broad surface of
neutral gas. An embedded source in the rim (Sugitani et al., 2002; Indebetouw
et al., 2007) may contribute to the 8 µm emission. At the sub-arcsecond
resolution of the 3.3 µm PAH image, the FUV-illuminated surfaces throughout
the pillar appear as composites of the many illuminated surfaces of sub-0.01
pc clumps. The wavy surfaces are reminiscent of the ridges and waves detected
by Berné et al. (2010) on the surface of Orion’s shell or by Hartigan et al.
(2020) along the edge of Carina’s Western Wall. The gas is clumpy/porous at
least down to this 0.01 pc length scale.

###  3.2 P1a: Threads and Cap

Figure 7: HCO+ (color) and CS (contour) line channel maps binned to 1 ${\rm
km~{}s^{-1}}$ and centered on the Pillars. The HCO+ images use an arcsinh
colorscale to make low-level emission more visible. The CS contours mark [0.7,
4.7, 8.7] K. The contour colors are chosen to increase contrast with the image
and have no further significance. The CS (left) and HCO+ (right) synthesized
beams are shown in the lower right corner of each map. The HCO+ (CS) 50% gain
contour is shown with a dashed (dotted) line. Figure 8: Position-velocity (PV)
diagram along each pillar using the [C II] and 12CO(J=1-0) line observations.
The left panel shows the F335N image with the three paths overlaid. The paths
are numbered and the labels are placed at the beginning of the path. The next
three panels show, from left to right, the PV diagrams of P1, P2, and P3 along
the numbered paths. In each PV diagram, the greyscale image is the CO PV
diagram at the native CO spatial resolution, while the contours show both CO
and [C II] at the [C II] resolution. Contour labels are main beam temperatures
$T_{\rm MB}$ in Kelvins. Emission at ${\rm V_{LSR}}\gtrsim 27~{}{\,{\rm
km~{}s^{-1}}}$ is from background features not directly related to the
Pillars. Velocities are VLSR. Figure 9: Same layout as in Figure 8, but with
PV diagrams along different paths. The first path runs across the Cap and
shows its peak velocity gradient. The second path crosses the two Threads and
shows the differences in [C II] and CO velocity structure. The third path
crosses the head of P2 and shows the spatial and kinematic offset of the [C
II] peak with respect to CO as well as the velocity gradient across the head
in both lines. Velocities are VLSR.

P1a appears composed of three morphologically and kinematically distinct
features: the “Cap” oriented across the top of Pillar 1 and two filamentary
structures (henceforth “Threads”; Figure 1) hanging down from the head of the
pillar in the NIR images in absorption (JWST F090W and F187N in Figure 1 and
the 1–2 µm images by Thompson et al. 2002) as well as in most molecular lines
in emission. The HCO+ and CS line channel maps in Figure 7 show that the
Threads are more or less parallel to each other and have similar lengths of
$\sim$0.3 pc. There are red-to-blue radial velocity gradients of order $-1$
${\rm km~{}s^{-1}}$ pc-1 along each Thread towards the pillar head (see
between 70–100″ and ${\rm V_{LSR}}<27~{}{\,{\rm km~{}s^{-1}}}$ in the
position-velocity diagram along P1 in Figure 8). The Eastern Thread is
spatially broader, $\sim$0.15 pc wide in CO emission just below the head of
the pillar, and is redshifted by $\sim$1 ${\rm km~{}s^{-1}}$ with respect to
the thinner ($\sim$0.07 pc) Western Thread. The Eastern Thread has a
transverse gradient in radial velocity, evident in the 12CO position-velocity
(PV) diagrams in Figure 9, redshifting to the west. [C II] emission towards
the Threads is centered on the Eastern Thread and appears uniform across that
section of the pillar rather than threaded like the molecular gas; this
difference is demonstrated in Figure 3 and is still evident when molecular
lines are convolved to the [C II] beam. A detailed analysis of the velocity
structure in [C II] and HCO+ emission, described in Appendix C, reveals a
kinematic detection of the Western Thread in [C II] spectra. Bonne et al.
(2023a) observes filamentary features, which they termed “legs”, extending
from the IC 63 nebula away from the stars with similar velocity gradients.

The peak brightness temperatures of molecular lines towards each Thread are
comparable. This implies one or both of 1) higher column density towards the
Western Thread, which in turn implies higher density assuming cylindrical
symmetry, or 2) higher temperature towards the Eastern Thread. The N2H+ line
emission towards P1a around ${\rm V_{LSR}}=23$–26 ${\rm km~{}s^{-1}}$ extends
slightly towards the northernmost part of the Western Thread, indicating
higher density.

The Cap is a lower velocity (${\rm V_{LSR}}=22$–24 ${\rm km~{}s^{-1}}$) [C II]
component associated with P1a which extends east from the head and lies on top
of the pillar like a cap. The Cap is observed in the lines of [C II], [O I],
and all but the highest critical density molecular lines. A dark NIR feature
indicating high extinction lies towards the Cap’s molecular line emission,
while the [C II] and [O I] are coincident with a bright NIR/illuminated-PAH
rim spatially shifted towards the exciting stars.

There is a steep gradient in radial velocity along the Cap, with velocity
increasing from east to west, in all lines in which the Cap is observed (see
the first PV diagram in Figure 9). Line widths broaden significantly towards
the center of the Cap where it meets the rest of P1a’s associated velocity
components. If we follow the Threads northwest into the pillar head, we find
that they spatially merge just as their velocity gradients lead their line
profiles to blend together in velocity and become indistinguishable. The line
profiles of both Threads and the Cap blend together just northwest of the
middle of P1a, where we detect the brightest molecular line and FIR continuum
emission. Detection of the N2H+ line towards this “Merge Point”, marked in
Figure 1, indicates cold, shielded gas. The [C II] emission on the eastern
side of the Eastern Thread has a particularly broad line profile with a lower
velocity wing. A similar pattern is observed in 12CO lines, and is present but
significantly diminished in HCO+ when compared to CS (see the spectra in
Appendix A).

###  3.3 P1b: Base and Horns

P1b, the lower half of P1, includes features from the Horns down to the
eastern half of the Shared Base. The complex of features is bright and
continuous in warm tracers like [C II], PAH emission, and FIR dust continuum
and discontinuous, particularly along the Shared Base, in molecular lines and
other dense gas tracers. The two Horns extend northwest from P1b in all
tracers with spatial resolution better than $\sim$10″. These $\sim$0.1 pc
diameter clumps are particularly pronounced in the molecular line channel
maps, indicating the presence of dense molecular gas, and the N2H+ and C18O
lines are faintly detected towards the Western Horn but not the Eastern Horn.
The Western Horn is brighter in longer wavelength ($>160~{}\mu$m) dust
continuum, while the Eastern Horn is brighter in NIR ionized gas tracers and
PDR tracers like [C II] and PAH features, consistent with a higher molecular
gas column density through the Western Horn. A bright rim of optical, NIR, and
PAH emission lies atop each Horn’s molecular emission in the direction of NGC
6611. We discuss the illumination structure of P1b further in Section 4.

We observe velocity gradients along each Horn: the top of the Western Horn is
more blueshifted than the gas below it, while the Eastern Horn is more
redshifted than the gas below it. Interpreting these Horns as pillar-like
structures along which gas is accelerated away from the stars (Pound, 1998),
the observed gradients imply that the Western Horn faces towards us and the
Eastern Horn faces away as we illustrate in Figure 10.

Below the Horns lies the Shared Base, a bright [C II]-emitting feature which
is distinguishable from the rest of P1 by its darker appearance in the
optical/NIR images (Figure 5). The high-resolution 3.3 and 8 µm PAH emission
maps show extended emission all the way across the Shared Base between P1 and
P2, indicating a broad illuminated surface. The Shared Base is the site of
broad line emission in all observed lines and a strong radial velocity
gradient ($\sim-3~{}{\,{\rm km~{}s^{-1}}}$ pc-1) smoothly connecting P1 to P2.
The [C II] and 12CO lines contain excess low-velocity emission between ${\rm
V_{LSR}}\approx 21$–22 ${\rm km~{}s^{-1}}$ compared to higher critical density
molecular lines like HCO+; we show example spectra in Appendix A.

###  3.4 P2

Continuum images from 3–500 µm show a cohesive column with a bright top. The
3.3 and 8 µm images resolve the head of P2 as a bright, $\sim$0.01 pc wide rim
atop a dark clump about 0.08 pc in diameter and reveal a second similarly
sized dark clump 0.4 pc below the head, about halfway down the pillar body.
This clump along P2’s body is outlined at its top by thin (0.006 pc at 3.3 µm)
bright rim and coincides with a re-brightening in the 70–500 µm images.

Line emission peaks between ${\rm V_{LSR}}=22$–23 ${\rm km~{}s^{-1}}$ towards
the pillar head; high density tracers like N2H+ peak closer to ${\rm
V_{LSR}}=23~{}{\,{\rm km~{}s^{-1}}}$ while warmer, lower density tracers like
[C II] and 12CO(J=1-0) are blueshifted by about 1 ${\rm km~{}s^{-1}}$ and all
12CO line profiles have low-velocity wings. This relative velocity shift of
warm tracers to dense tracers, similar to low-velocity wings observed towards
P1b and to lesser degree towards P3, persists throughout the northern half of
the Pillar, but the pattern disappears as line widths increase towards the
southern half. We discuss this behavior in Section 6.1.

Velocity gradients are detected in most lines both along and across the pillar
body (Figures 8 and 9). Gradients across P1 and P2 were also observed by Sofue
(2020). Line velocities are highest towards the northeastern side of the
pillar head, and decrease both to the west across the head and to the south
along the body. [C II], and [O I] where it is observed towards the head, trace
a coherent column which gets almost monotonically brighter towards the base of
the pillar in integrated intensity, with only a slight local brightness maxima
at the head where the PAH tracers brighten. Molecular lines trace more
substructure along P2 than [C II]. HCO+, HCN, CS and the CO lines trace an
elongated clump towards the head, the clump towards the middle of the body,
and two filametary tails below the second clump. The C18O and N2H+ lines are
only detected towards the head and the mid-body clump, where all molecular
line profiles are broadest. Between these two features, molecular line
emission is dim and lines are narrower (see for example the PV diagram along
P2 in Figure 8).

The NIR and optical (Hester et al., 1996) images feature a bright “wisp” about
midway along the pillar, just above the location of the clump. The wisp seems
to originate from P3 and cross the body of P2, as we see some continuous edges
in the optical images both towards and off P2. Emission from the HCO+, CS, and
13CO lines is dimmer where the wisp overlaps with P2 between the head and the
clump. The wisp, where it passes over P2, is therefore coincident with a
region of optically thinner lines of sight through P2. The southern edge of
the wisp is very close (0.05 pc or 6″, smaller than our molecular line beams)
to the dark, dense clump mid-way down the pillar body.

###  3.5 P3

Pillar 3, the smallest of the three main Pillars, presents in our molecular
line observations as a $\sim$0.38 pc long pillar oriented roughly parallel to
P1 and P2 with two tails which extend $\sim$0.34 pc in either direction at
$\sim$100∘ angles from the body in a “wishbone” shape also seen in optical
images. The head of the pillar is brightest at ${\rm V_{LSR}}\approx
21.2~{}{\,{\rm km~{}s^{-1}}}$ in all molecular tracers. Using the line profile
modeling described in Appendix C, we find that [C II] is blueshifted relative
to the molecular lines by $\lesssim 0.3~{}{\,{\rm km~{}s^{-1}}}$ through the
head. C18O and N2H+ are only observed towards the pillar head. The wishbone
tails are reminiscent of the “ears” of cometary globules simulated by Lefloch
& Lazareff (1994) in 2 dimensions, which they discuss in Sections 5.2 and
5.3.1 of their paper.

###  3.6 Ridge

The Ridge and P4 (Section 3.7) are detected in most molecular lines and are
each associated with a compact clump of N2H+ line emission indicating enhanced
gas density. The Ridge is outside of the half-power beam of the CO(J=1-0)
observations, and P4 was not covered by the CO(J=6-5) observations and lies on
the edge of the HCO+, HCN, CS, and N2H+ half-power beam.

The Ridge, lying around 0.5 pc southeast of the main Pillars’ tails and
appearing in the same ${\rm V_{LSR}}\approx 25$ ${\rm km~{}s^{-1}}$ channel
maps as P1 (Figure 7), is a 0.8 pc long and 0.08 pc thick bar of molecular and
[C II] line emission oriented perpendicular to the direction of illumination.
The Ridge spans nearly the same width as the three-pillar system. The
integrated FIR dust emission from the Ridge (proxy for FUV radiation field;
Section 4) is shifted towards NGC 6611 by $\sim$0.03 pc with respect to dust
column density (Section 5), and at higher resolution, the 3.3 and 8 µm PAH
emission arise from the edge of the Ridge facing the cluster and decrease in
brightness towards sites of enhanced molecular line emission along the Ridge.
We do not detect a spatial shift of [C II] line emission with respect to
molecular line emission.

Radial velocity increases from east to west along the Ridge in all lines in
which it is observed. Below the Ridge, we observe diffuse PAH and dust
emission between 3–500 µm as well as diffuse [C II] emission between ${\rm
V_{LSR}}=24$–27 ${\rm km~{}s^{-1}}$ (same velocity interval as the Ridge)
clearly bounded to the northwest by the Ridge and to the northeast by another
ridge of gas perpendicular to the Ridge. This diffuse emission appears bounded
to the west/southwest by a curved stream of gas, represented by the dashed
line in Figure 1, but this stream of gas is separated by a few ${\rm
km~{}s^{-1}}$ from the Ridge in [C II], CO(J=3-2), HCN, and HCO+ line
velocity. To the south, the diffuse emission continues into the bright,
$\sim$10 pc scale feature associated with the edge of the H II region. We do
not detect this diffuse emission in the molecular lines, but we do detect the
Ridge, the perpendicular ridge to the northeast, and, faintly, the
west/southwest stream. The Ridge and northeast ridge are comparably bright in
[C II], while the northeast ridge is somewhat dimmer in 12CO(J=3-2), much
dimmer in HCO+ and HCN, and not detected in CS.

Molecular gas column density is higher along the boundaries between the
diffuse neutral gas and the ionized gas, particularly the Ridge, and all the
neutral gas between the Ridge and northeast ridge is well illuminated. The
west/southwest stream is closer in velocity to P2 and P3, while the Ridge is
more kinematically similar to P1 as shown by the velocity RGB image in Figure
5. Figures 1 and 5 show optical and NIR counterparts to the Ridge, northeast
ridge, and west/southwest stream.

###  3.7 P4

P4, a more triangular feature than the three main Pillars, has bright eastern
and western edges in the 3.3 and 8 µm PAH maps. NIR continuum, like that in
the central panel of Figure 1, shows extinction in all bands through the head
and body. From the [C II] channel maps, we detect a radial velocity gradient
along the edge of P4, from its southwestern corner (${\rm V_{LSR}}\sim
21~{}{\,{\rm km~{}s^{-1}}}$) up its western side to its point ($\sim
23~{}{\,{\rm km~{}s^{-1}}}$) and back down its northeastern corner ($\sim
25~{}{\,{\rm km~{}s^{-1}}}$). [C II] line emission outlines the (upper two)
triangular edges of P4 and is around 30% dimmer towards the middle of the
pillar. The blueshifted western edge corresponds with the dark ridge in the
NIR images. We do not resolve much spatial structure in the CO(J=3-2)
observations but do detect some variability in radial velocity along the
feature. P4 lies towards the edge of the H II region, which is a brighter [C
II] source with a complex line profile that is not well separated from P4’s
line emission.

###  3.8 Summary and Geometry

The overall picture that emerges from this wealth of data is schematically
represented in Figure 10. The PAH emission traces surface structures
illuminated by the stellar cluster and can be used to derive the geometry of
the Pillars and their orientation relative to the stellar cluster. The
Pa-$\alpha$ emission originates mostly from diffuse ionized gas filling the H
II region surrounding the Pillars. The relative strength of this line helps in
placing structures such as P4 and the Ridge more on the near side or far side
of the ionized cavity. Peak [C II] and molecular line velocity gradients along
each pillar trace bulk motions of the gas pushed away from the stars and can
indicate whether pillars are inclined towards or away from the observer
(Pound, 1998; McLeod et al., 2015). This analysis places P1b, P2, and P3 in
the foreground and, while P1b and P2 point toward the illuminating stars, P3
is backlit. P1a is in the background and points toward the illuminating stars
as well. P1b and P2 are connected by the Shared Base, while a connection
between P1a and P1b is implied by the kinematically continuous [C II] emission
between them and their on-sky projected alignment. Warm gas tracers are
blueshifted by $\sim$0.5 ${\rm km~{}s^{-1}}$ towards P1b and the head of P2.
Since these pillars are both facing away from us, we interpret this to mean
that the less dense outer layers are accelerated down the pillars more quickly
than the dense interior layers; this is discussed in detail in Section 6.1. We
note that we determine the sign of the pillar inclination (towards or away
from the observer), but not its magnitude, from the kinematic analysis; the
particular inclinations depicted in Figure 10 are but one possible
configuration. Sofue (2020) calculates the absolute values of the inclinations
to be 47$\arcdeg$, 40$\arcdeg$, and 40$\arcdeg$ for P1a, P2, and P3,
respectively; each may be towards or away from the observer.

The Threads’ radial velocity gradients, redshifting away from the head, are
consistent with material being accelerated away from the dense head along a
pillar which is on the far side of the illuminating stars and points towards
both the stars and the observer, consistent with the LOS position of P1a
suggested by Pound (1998) and McLeod et al. (2015). The Eastern Thread is
redshifted with respect to the Western Thread; if material flows down both
Threads with a similar velocity, then the Eastern Thread may be more inclined,
as drawn in Figure 10, so that it has a greater projected velocity than the
Western Thread. The Cap’s radial velocity gradient blueshifts away from the
head, which would imply that it lies between the stars and the observer and
faces away from us. The velocity gradient may be due to unknown kinematic
interactions in the pillar head, or it is possible that the head is extended
along the line of sight so that parts of it are on either side of the cluster
with respect to the observer. This second explanation is not so unreasonable
if we expect that P1’s head is almost right below the cluster with very little
line-of-sight offset; in this case, the line-of-sight separation between the
nearest and farthest parts of the head, on opposite sides of the cluster with
respect to us, wouldn’t be that large. We elect to place P1a directly below
and only slightly behind the stars in our geometric model in Figure 10, in
accordance with the second explanation for this velocity gradient.

The Shared Base is depicted in Figure 10 extending a significant distance
between P1 and P2. Our analysis of column density towards the Shared Base
versus average number density elsewhere implies a LOS width $\gtrsim$0.5 pc
(see Sections 4.1.2 and 5.1), and our analysis of the FUV radiation field
places the heads of P1a and P2 $\sim$1 pc from each other along the LOS (see
Section 4.2). Both of these analyses depend on the geometry and optical depth
of the Pillars, and the latter is sensitive to the LOS geometry of the cluster
members and extinction of FUV radiation between the cluster and the Pillars,
though relative distances (such as between P1a and P2) via the latter method
are not affected by uniform extinction between the cluster and all features.
Figure 10 represents our analysis-based educated guess of LOS widths and
separations of the Shared Base and other features, but these uncertainties
prevent us from making precise estimates.

The positions of P4 and the Ridge are more uncertain than those of the three
main Pillars. We place P4 closer to the observer along the LOS based on its
relative darkness in Pa-$\alpha$, and the Ridge farther from the observer
based on its relative brightness in Pa-$\alpha$ (Figure 1). FIR emission
towards both of these features, but particularly P4, is brighter than we
expect given the projected distances of the selected stars, so the assumptions
about the cluster which we use to estimate incident FUV radiation field for
the three primary Pillars are not appropriate for features farther away from
the cluster. Additionally, some of P4’s FIR emission might originate from 1)
an embedded young stellar object (YSO) and/or 2) internal heating by that YSO,
rather than solely from reprocessed FUV radiation from the main cluster. We
discuss this in more detail in Section 4.2.

Figure 10: Schematic diagram of the pillar system as viewed from the
southwest, so that the observer’s line of sight runs horizontally across the
figure. The green lines mark molecular gas structures and the yellow highlight
marks atomic gas structure, which forms an envelope around the molecular gas
(Section 4.1.1). The gold rays originating from the top of the figure mark the
direction of illumination, with the solid central line marking the
“perpendicular” ray and the assumed line-of-sight position of the stars. P1a
and P1b are both collections of features. The Shared Base and P1b overlap in
definition, and the distinction between them is not particularly relevant
because they are both approximate labels for (sets of) observed features. We
estimate that the heads of P1a and P2 are separated by $\sim$1 pc along the
line of sight; see Section 3.8 for more detail.

##  4 Photodissociation Regions

###  4.1 Morphology and Geometry of the Major PDRs

Any bright [C II]-emitting feature is likely the site of a PDR, and our [C II]
observations show that M16 is rich in PDR emission. The vast majority of the
[C II] emission is indeed from atomic gas, with a line width of $\sim$1 ${\rm
km~{}s^{-1}}$ rather than from ionized gas which would have a turbulence-
dominated line width closer to 10 ${\rm km~{}s^{-1}}$ (Cuadrado et al., 2019,
and see Appendix B for discussion about [C II] from the H II phase).
Throughout P1 and the rest of the pillar system, sub-arcsecond resolution
optical and near-IR maps trace bright rims which are $\sim$5–20 times brighter
than their surroundings. We interpret these areas, which include the Cap,
Eastern Horn, and Shared Base in P1, the bright rim atop the head of P2, and
parts of the head of P3 and the Ridge, to be limb-brightened edge-on PDRs. We
interpret the dimmer, “average brightness” emission to originate from poorly
illuminated PDRs or PDRs viewed face-on or through dust extinction.

The remainder of this section highlights some of the bright PDRs throughout
the Pillars with a focus on their orientation with respect to the observer and
implications about the overall geometry of the system. In order to discuss the
geometry of some of these structures, we must consider the atomic and
molecular gas column densities through them. Column densities are discussed in
depth in Section 5, but the results of that section inform the following
discussions.

####  4.1.1 P1a

Bright PDRs are associated with three distinct sets of features along the body
of P1: the Cap, the Horns, and the Shared Base. The head of P1 features a
prominent PDR across the top of the Cap, observed as a bright rim in 8 µm and
[C II], and [O I] sitting atop the molecular emission (Figure 3), which we are
viewing edge-on (Levenson et al., 2000). Below the Cap, we see relatively
bright PAH emission towards a sort of “shoulder” atop the Western Thread, but
only moderate PAH emission between the Cap and the Horns indicating that there
are few edge-on illuminated surfaces towards the Threads.

Most of the [C II] and molecular spectra towards the head of P1a are not fit
well by a single Gaussian component. Since we observe several morphologically
distinct components in the channel maps presented in Section 3.2, we determine
that P1a contains multiple components which are separated by less than their
individual linewidths. In a detailed kinematic analysis of P1a described in
Appendix C, we find that the HCO+ spectra towards P1a are fit well with 3
components corresponding to the Cap, Eastern Thread, and Western Thread, while
the [C II] spectra towards P1a are fit well with 2 components corresponding to
the Cap and a combined-Thread component. The [C II] spectra towards the
Threads are dominated by the Eastern Thread component, which lies at a higher
velocity than the Western Thread component based on their molecular line
velocities, but we detect a weak signature of the Western Thread component in
the [C II] spectra. As the [C II] emission does not exhibit the threaded
morphology found in the molecular emission, even at matched spatial
resolution, we determine that the [C II] likely originates from a more
uniform, extended envelope of atomic gas surrounding the dense molecular gas
features.

Unlike the P1a spectra discussed above, the [C II] and molecular spectra
towards the Merge Point in P1a, marked in Figure 1, are fit well with a single
component. Based on our analysis in Appendix C, we suggest that the Cap and
two Threads are physically joined together as one component towards this
position. N2H+ line emission detected towards the Merge Point indicates cold,
shielded molecular gas buried deep within the cloud. Our proposed geometry in
Figure 10 envisions the Cap as a compressed rim or globule-like structure and
the two Threads as legs trailing down from that cloud.

####  4.1.2 P1b

The Horns, the warm gas below them, and the Shared Base are all capped with
bright rims of PAH emission and form a terraced arrangement of edge-on PDRs.
These rims lie atop sites of molecular emission which are spread out in
velocity. Molecular line emission towards the Horns lies at ${\rm
V_{LSR}}=23$–25 ${\rm km~{}s^{-1}}$. The bases of the Horns are connected by
warm, low column density 12CO-bright gas between ${\rm V_{LSR}}=23$–24 ${\rm
km~{}s^{-1}}$ (White et al., 1999). The Shared Base lies at ${\rm
V_{LSR}}=22$–23 ${\rm km~{}s^{-1}}$. The Shared Base and the Horns are
spatially and dynamically connected but are separated by about 0.13 pc and 1
${\rm km~{}s^{-1}}$ and appear as distinct features in molecular line channel
maps. These observations are consistent with the Horns and Shared Base being
distinct, but connected, sites of strong PDR activity. P1b and the Shared Base
host multiple clumps of illuminated gas and behave in some ways like pillar
heads.

Table 3 in Section 5 lists a high atomic column density $N({\rm H})$ derived
from C+ towards the Shared Base compared to what we expect based on the PDR
modeling described in Section 4.3. Either the atomic gas density towards the
Shared Base is $\sim$2 times higher than other regions or there is a long
line-of-sight length through this gas structure. We prefer the second
explanation since we do not expect P1b to be the site of enhanced atomic gas
density, particularly since there is little molecular gas. If we expect the
density to be closer to the median atomic gas density in Table 3, then the
Shared Base may extend $\gtrsim$0.5 pc along the line of sight like a “valley”
between P1b and P2 as pictured in Figure 10. This can explain the strong
observed gradient in peak velocity along the Shared Base, as the projected
gradient would be larger than the physical gradient. The terraced arrangement
of illuminated surfaces along P1b as well as the position of P1a below NGC
6611 while P1b and P2 are closer to the observer (Pound, 1998; McLeod et al.,
2015) suggest that P1b continues to extend away from the observer along the
LOS above the Shared Base, towards the Horns.

The emission characteristics of the Eastern and Western Horns described in
Section 3 indicate that the Eastern Horn is associated with more illuminated
surface, either in total or facing us. Some, but not all, of the Eastern
Horn’s brightness compared to the Western Horn can be explained by two
geometric phenomena. First, the Western Horn’s illuminated surface must be
mostly on the far side of the feature so that its $<$8 µm emission is
extincted by dust within the Horn. Second, the Eastern Horn is superimposed on
the body of P1, which can still be viewed through the Horn at longer
wavelengths (dust continuum $\geq$70 µm and the [C II] line) and must be
responsible for the Eastern Horn’s brightness at those wavelengths. This line-
of-sight relationship is easiest to see in the optical images where it is
clear that there is extra pillar emission surrounding the Eastern Horn. The
observed radial velocity gradients (Section 3) along the Horns suggest that
the Western Horn faces towards us and the Eastern Horn faces away.

####  4.1.3 P2

All PDR tracers (PAH, 70–160 µm, [C II]) are dimmer towards P2 and P3 than
towards the PDR-heavy P1. The brightest [C II] emission along P2 lies around
${\rm V_{LSR}}=20$–22 ${\rm km~{}s^{-1}}$ and is associated with the Shared
Base. The emission further up the body of P2 is roughly 60% as bright as the
emission towards its base and remains roughly constant in brightness even
towards the head of P2. We attribute these observations to P2 hosting a
smaller illuminated PDR surface than the various locations along P1. Towards
P2, we observe PDR tracers through a limited column density and thus they
appear less bright than towards P1. Bright PAH emission is observed along the
edge of the P2 head, which is not as flat in projection as the Cap in P1a and
may present 1) less surface area in total and 2) a larger fraction of surface
area illuminated at higher inclination (greater angle from the normal). Below
the head of P2, we see a rim of enhanced PAH emission atop the dark clump mid-
way down the body and atop several cometary clumps towards the base. P1 and
P2, despite a clear difference in the brightness and abundance of PDRs along
their bodies, both contain numerous small illuminated surfaces between their
heads and tails.

####  4.1.4 P3

P3 is even dimmer in [C II] than P2 (see spectra in Appendix A) and is not as
well resolved in [C II] or 70 µm since its projected width ($\sim$14″, 0.1 pc)
is comparable to their beams. At 3.3 µm, the head of P3 is capped by a thin
($\approx$0.008 pc) rim of PAH emission. The two tails of P3 which give it its
wishbone shape are brighter along their cluster-facing edges in 3.3 and 8 µm
emission, indicating that the tails host PDRs. As P3 and its tails are small
compared to the beams of most of our observations, we will not be able to
study them in as much detail as P1 and P2.

####  4.1.5 Ridge

The fact that the Ridge is approximately perpendicular to the direction of
illumination and spans the width of the three-pillar system lying above it
suggests that it may have formed in the shadow of the Pillars. It is at least
somewhat illuminated by the cluster at present, as it hosts an extended PDR
surface along its 1 pc long surface facing NGC 6611. The derived C+ column
density towards the Ridge is higher than expected given its projected width
($\approx$0.1 pc), similar to what we described for the Shared Base in Section
4.1.2, so it must have a LOS size of $\gtrsim$0.3 pc. The rim of PAH emission
is brighter and thicker where molecular line emission is weak along the Ridge,
which could be due either to extinction through the dense clumps if the Ridge
faces slightly away or gas density variations along the Ridge altering the
physical PDR width if the Ridge is viewed nearly edge-on.

####  4.1.6 P4

The two sides of the angular P4 both host edge-on PDRs at their surfaces. We
assume the low 1.87 µm Pa-$\alpha$ emission towards P4 (green in the central
panel of Figure 1) is due to a low foreground column density of diffuse
ionized gas, indicating that P4 is relatively close to the observer within the
pillar system. The low intensity PAH emission towards the face of P4 may then
originate from the illuminated far side and be extincted by dust inside the
structure. We cannot determine the orientation of P4 based on kinematics or
extinction, but we orient it in Figure 10 pointing away from the observer
towards the main cluster. One must consider the possibility that it faces
towards the observer similar to P3.

The source driving the Herbig-Haro object HH 216 lies at the tip of P4
(Andersen et al., 2004; Indebetouw et al., 2007; Flagey et al., 2020).
Indebetouw et al. (2007) suggest that it is one of two Class I YSOs identified
towards the tip of P4, while Flagey et al. (2020) identify a third nearby
point source which may be responsible. We discuss the possible effect of this
source on the FUV radiation field at P4 at the end of Section 4.2.

###  4.2 Sources of Illumination

The brightest members of NGC 6611 are $\sim 2.5$ pc from the Cap in P1a and
$\sim 5$ pc from P4, meaning that illumination varies by a factor of $\sim 4$
along the pillar system. We use the catalog of NGC 6611 published by
Hillenbrand et al. (1993) to estimate the FUV radiation field $G_{0}$
throughout the pillar system. The process is described in Tiwari et al. (2021)
and the software is publicly available as scoby555The code is archived at
http://hdl.handle.net/1903/30441; scoby is also developed on GitHub. We link
the spectral types of early-type NGC 6611 members published by Hillenbrand et
al. (1993) to $T_{\rm eff}$ and ${\rm log}~{}g$ values using the tables of
Martins et al. (2005), and then use the $T_{\rm eff}$ and ${\rm log}~{}g$
values to select stellar models from the PoWR suite (Sander et al., 2015;
Hainich et al., 2019). The PoWR stellar models provide theoretical spectra,
which we combine with the catalog coordinates to estimate the FUV (6–13.6 eV)
intensity at a given location from a single star assuming the projected
distance and no intermediate extinction. We estimate the total FUV intensity
at the given location by summing across the contribution from all, or a subset
of, stars in the catalog. There are a handful of early-type stars within
$\sim$1 pc of the Pillars in projection which cause our $G_{0}$ estimate to
exceed 5000 Habing units towards the Pillars. This is unreasonable given the
FIR-based estimates described later in this section, so these stars must lie
$\gtrsim$1 pc away along the line of sight. We restrict our sample to stars
within 2 pc in projection of the approximate cluster core
$(\alpha,\,\delta)_{J2000}=(18\mbox{${}^{\mathrm{h}}$}18\mbox{${}^{\mathrm{m}}$}35\fs
9543,\ -13\arcdeg 45\arcmin 20.364\arcsec)$, which includes 30 stars with
types B2.5 and earlier. We additionally restrict stars to have ${\rm
log}_{10}(L_{\rm FUV}/L_{\odot})>4.5$ since the largest cluster members should
dominate the feedback; 8 stars remain after the filtering. We find that
including all 30 stars increases $G_{0}$ by $\sim$50% near the Pillars. The
FUV radiation field calculated from our final list of 8 stars (# 161, 166,
175, 197, 205, 210, 222, 246 in the catalog of Hillenbrand et al. 1993) is
$G_{0}\sim 2500$ Habing units towards the heads of P1a and P2, $\sim 1700$
towards P3, $\sim 1500$ Habing units towards the Shared Base, and $<1000$
Habing units below the Ridge.

We can also estimate the FUV radiation field based on FIR thermal emission
from dust, assuming that dust is re-radiating absorbed FUV radiation (Wolfire
et al., 2022). The FUV radiation field values are derived from the FIR dust
emission at 70 and 160 µm using the methodology described in Schneider et al.
2016 and corrected for a background sampled from the regions described in
Appendix B. Through this method, we estimate $G_{0}\sim 2000$, 700, and 300
Habing units towards the heads of P1a, P2, and P3, respectively; and 800–1200
Habing units towards the Shared Base, the Ridge, and P4.

The stellar estimates of $G_{0}$ should be strict upper limits as they assume
that all stars and features lie at their projected distances on the plane of
the sky and that all starlight reaches every feature without extinction. The
discrepancy between the higher stellar estimate and lower FIR estimate of
$G_{0}$ towards P2 and P3 could be a result of 1) line-of-sight separation on
the order of a few parsecs, consistent with their positions on the near side
of the cluster with respect to the observer, 2) intermediate extinction by
gas/dust within the cavity, and/or 3) beam dilution by the $\sim 14$″ beam at
160 µm, as the PDR surfaces at the heads of these two smaller pillars cover a
smaller solid angle than the PDR towards the Cap in P1. Allowing for some
line-of-sight spread among the cluster members, the FIR-based estimates of
$G_{0}$ towards the Cap, Shared Base, Ridge, and P4 are fairly consistent with
the stellar estimates. These features must lie close to their projected
distances from the stars; this is most significant for the Cap, which lies
fairly close in projection and so must lie approximately directly below the
brightest cluster members.

The Pillars’ FUV illumination is dominated by the massive cluster members near
the cluster core. The handful of early-type stars which are close in
projection (within $\sim$1 pc) to the Pillars are not significant contributors
to the FUV field at the Pillars and thus must be separated from the Pillars by
$\gtrsim$1 pc along the line of sight, though they may still be important
contributors to the FUV field below the Ridge and to off-axis illumination.
These near-Pillar stars may shine on parts of the pillar bodies that would
otherwise be shadowed by the pillar heads with respect to the cluster core;
that said, off-axis illumination is probably dominated by the diffuse extreme
ultraviolet (EUV, $h\nu>13.6$ eV) field created by recombinations directly to
the ground state.

P3 appears not to face the bright cluster core, as explained in Section 1. It
is possible that it faces a different star than the rest of the pillars, one
which lies on the outskirts of the cluster. P3 is the nearest of the Pillars
to the observer (McLeod et al., 2015, and see schematic in our Figure 10), so
a star displaced closer to the observer along the LOS could have a significant
impact on the evolution of P3 while having little impact on P1 or P2. We do
not have sufficient information to identify a candidate star, but suggest that
it could be among stars # 166, 197, 205, and 210 (O8.5 V, O7 V((f)), O5
V((f*)), and B1 III, respectively), which lie $\sim$2 pc away in projection
from P3 in the direction it points. P1 and P2 are well aligned with the
expected morphology and kinematics of the radiative interaction of the NGC
6611 cluster with the molecular shell; we consider it unlikely that P3’s
orientation away from the cluster is evidence that P3 was sculpted by a
different phenomenon than the other Pillars, such as ablative Rayleigh-Taylor
instability (Mizuta et al., 2006), given its similarity in morphology and
close proximity. Since we cannot constrain the orientation of P4, we must
consider the possibility that it faces us like P3, in which case this same
discussion would apply.

####  4.2.1 Line-of-Sight Pillar Geometry from $G_{0}$

We can roughly estimate the line-of-sight geometry of the Pillars using our
two $G_{0}$ estimates, assuming that the catalog-based $G_{\rm 0,\,star}$ is
the emitted FUV radiation field and the FIR-based $G_{\rm 0,\,dust}$ is the
apparent FUV radiation field seen by the features. The stars are assumed to be
at their projected distances, lying on the plane of the sky, and all optical
paths between the cluster and the features are assumed to be free of
extinction. This model attributes all differences between $G_{\rm 0,\,star}$
and $G_{\rm 0,\,dust}$ to the line-of-sight displacement of each feature from
the plane on which the cluster stars lie. To simplify the calculation, we
further assume that all the stellar radiation originates from the center of
the cluster using the same coordinate as earlier. The catalog-based estimate
assumes that all features are at their projected distances, so $G_{\rm
0,\,star}=L/4\pi r^{2}$, where $r$ is the projected distance from a given
position to the cluster center and $L$ is the FUV luminosity of the cluster.
The FIR-based estimate measures the apparent radiation field, so $G_{\rm
0,\,dust}=L/4\pi(r^{2}+z^{2})$, where $z$ is the displacement of the feature
from the cluster plane. The ratio of the two $G_{0}$ estimates can be solved
for the absolute value of $z$.

$\frac{G_{\rm 0,\,star}}{G_{\rm
0,\,dust}}=\frac{r^{2}+z^{2}}{r^{2}}=1+(z/r)^{2}\ \rightarrow\
|z|=r\Bigg{(}\frac{G_{\rm 0,\,star}}{G_{\rm 0,\,dust}}-1\Bigg{)}^{1/2}$ (1)

The heads of P1a and P2 are estimated to lie $\sim$0–1 pc and $\sim$2–3 pc,
respectively, from the plane of the cluster; higher values are obtained when
the background subtraction is applied to $G_{\rm 0,\,dust}$. The head of P3 is
close to the size of the 160 µm beam, so beam dilution may slightly reduce
$G_{\rm 0,\,dust}$ and inflate separation; we estimate a separation of
$\sim$2–5 pc from the plane. The absolute separations for each feature are
sensitive to the background subtraction, which is in turn sensitive to the 70
µm opacity through the head, and also the underlying assumption of extinction-
free optical paths. The relative separation between features is less sensitive
to these. Given the assumptions and uncertainties involved, we conclude that
P1a/P2 and P2/P3 are each separated on the order of $\sim$1 pc along the line
of sight, so that the three primary Pillars span a few pc along the LOS.

We compare to the distances calculated by Sofue (2020) using peak radio
continuum brightness. LOS separations from the plane of the cluster $z=1.8$,
2.5, and 2.8 pc for P1a, P2, and P3 are calculated as $z=D{\rm cos}(i)$ from
the true distances $D=2.6$, 3.2, and 3.6 pc and inclinations $i=47\arcdeg$,
40$\arcdeg$, and 40$\arcdeg$ presented by Sofue (2020). In both their and our
estimates, P1a is nearest and P3 is farthest from the plane of the cluster.
Our estimate, from both $G_{0}$ and kinematic clues (Section 3.8), positions
P1a closer to the plane of the cluster than the estimate by Sofue (2020).

Line of sight displacement $z$ can only be calculated when $G_{\rm
0,\,star}\geq G_{\rm 0,\,dust}$, which is the case for the three primary
Pillars but not towards the Ridge or P4. Stars $>2$ pc from the cluster core
in projection may contribute more significantly to the radiation field at
those two features since they are farther from the cluster core. The LOS
displacements of the Ridge and P4 therefore cannot be estimated using this
method. The Shared Base is extended along the LOS (Section 4.1.2), so $G_{\rm
0,\,dust}$ is likely overestimated there and LOS distance cannot be
determined.

At the tip of P4 lie two Class I YSOs, HH-N and HH-S, identified by Indebetouw
et al. (2007) in NIR and mid-IR images, and a third source of unknown type
$\sim$3″ from HH-N identified by Flagey et al. (2020) in the NIR. Any one of
these sources may drive HH 216. A point source appears at the tip of P4 in
both the 70 and 160 µm observations near the location of HH-N and the NIR
point source. The source emission is $\sim$50% brighter than the diffuse P4
emission in both bands and likely affects the $G_{\rm 0,\,dust}$ estimate
towards that position. It is unclear whether radiation from the source heats
or otherwise affects P4. Indebetouw et al. (2007) noted that the YSOs appear
extincted by tens of $A_{V}$, but the NIR source’s luminosity and extinction
are not known. If P4 is heated internally, the derived $G_{\rm 0,\,dust}$ is
unsuitable for the LOS separation estimate.

###  4.3 Modeling PDRs Towards the Pillars

We model observations towards several PDRs associated with the Pillars to
measure densities and FUV radiation fields at those locations (see Wolfire et
al. 2022 for an overview of PDR modeling). We use the Wolfire-Kaufman 2020
models available in the PDR Toolbox which we access via the
pdrtpy666https://dustem.astro.umd.edu software (Kaufman et al., 2006; Pound &
Wolfire, 2008, 2023). This particular model set assumes a plane-parallel,
face-on PDR geometry most appropriate for PDRs at the surface of non-clumpy
clouds. The observed intensities of [C II], [O I], 12CO(J=1-0), 12CO(J=3-2),
and 12CO(J=6-5) at selected locations towards the PDRs are integrated over the
relevant velocity intervals for each PDR and convolved to the [C II]
resolution. We do not convolve 12CO(J=3-2), which has a slightly larger beam.
The higher resolution of the other observations is more valuable than keeping
all the beams matched.

The on-sky widths of the modeled features are close to the [C II] beam width,
so we assume that 1) the beam filling factor is 1 for all observations except
12CO(J=3-2), 2) there is zero emission outside the features and 3) the
12CO(J=3-2) beam filling factor is the ratio of the beam areas
$\Omega_{{{\rm[C\ II]}}}/\Omega_{\rm CO}$. We scale the 12CO(J=3-2)
intensities up by the ratio $\Omega_{\rm CO}/\Omega_{{{\rm[C\ II]}}}$ to
account for beam dilution.

For pairs of these integrated intensities, the PDR Toolbox takes ratios and
plots them on a grid of density $n$ and FUV radiation field $G_{0}$ using
precalculated models. This procedure cancels out the beam filling factor to
first order. The resulting diagrams, shown in Figure 11, are called overlay
plots and their usage is described in detail in Tiwari et al. (2022) and Pound
& Wolfire (2023). The ratios appear as lines curving across the grid. In
principle, all lines should intersect at a density and radiation field
consistent with the observations. In practice, our overlay plots generally do
not converge at one specific location but rather a region of the $n$ vs.
$G_{0}$ grid around $n\approx 2\times 10^{4}$ cm-3 and a few $10^{2}$ Habing
units, which we take to be the physical conditions consistent with our
observations under the Wolfire-Kaufman 2020 PDR models. The results of the PDR
Toolbox fitting routine for the four sample locations from Figure 11 are
listed in Table 2.

Density is fairly well bounded under this model, with our observations fitting
between $\sim$1–4 $\times 10^{4}$ cm-3. We estimate this uncertainty using a
contour of a few times the minimum $\chi^{2}$ and confirm that this region of
parameter space doesn’t change much between lines of sight which should have
similar density and illumination. This density estimate is consistent with the
estimates made by Pound (1998) and Levenson et al. (2000) as well as broadly
consistent with our estimate based on the [C II] column density and assumed
line-of-sight geometry (Section 5.1.2; Table 3).

Radiation field is much more loosely bound, as we can only place an upper
limit of $G_{0}\sim 1000$ Habing units and no realistic lower limit under this
model. The PDR models used here assume face-on geometry, which may produce
slightly inaccurate results, so we rely on independent estimates of density
and radiation field to gauge how accurate the PDR model-based estimates might
be.

We can further bound the radiation field seen by the Pillars using our
independent estimates of $G_{0}$ from stellar catalogs and FIR dust emission
which we described in Section 4.2 The stellar method, as previously discussed,
yields estimates of $G_{0}\sim 1$–2 $\times 10^{3}$ Habing units in the
neighborhood of the Pillars. These estimates mark the upper end of the PDR
model estimate of radiation field, so we assume that the radiation field seen
by the Pillars is somewhere between 500–2000 Habing units depending on local
geometry.

We are unable to make substantial claims about spatial variation in the
density and radiation field based on the PDR models, as the uncertainties in
our modeling techniques are greater than the local variation. We find that the
densities and radiation fields are generally self-consistent across all
modeled regions towards the Pillars.

Figure 11: (Top row) Overlay plots made using pdrtpy which show observed line
ratios towards 4 locations which are labeled in the top left corner of each
plot. The black cross shows the automatically fitted solution and associated
uncertainty; these values are listed in Table 2. (Bottom row) The reduced
$\chi^{2}$ associated with the automatic fit, with overlaid labeled contours.
Radiation field values estimated via stellar catalog and FIR emission are
plotted as horizontal dashed lines in each figure. Table 2: Fitted $G_{0}$ and
$n$ to observed line ratios using pdrtpy.

Location | $G_{0}$ | $n$ | $\chi^{2}/{\rm dof}$ | dof
---|---|---|---|---
Name | (Habing) | (cm-3) |  |
NE-thread | 150 | $2.6\times 10^{4}$ | 150 | 5
W-Horn | 460 | $2.7\times 10^{4}$ | 14 | 5
P2 | 300 | $3.1\times 10^{4}$ | 58 | 5
P3 | 89 | $2.3\times 10^{4}$ | 6.3 | 4

Note. — The $G_{0}$ and $n$ values are those fitted to the observed line
ratios overlayed on the top row of Figure 11. These correspond to the minimum
$\chi^{2}$ location, in color on the bottom row of Figure 11, which is
indicated with a black cross in all panels of that figure.

##  5 Mass and Physical Conditions

###  5.1 Molecular and Atomic Hydrogen Column Densities

Molecular and atomic gas column densities are estimated using the 13CO(J=1-0)
and [C II] lines, respectively, as well as with dust emission. Column
densities are summed pixel-by-pixel to estimate the total pillar masses, and
line-of-sight distances are assumed to estimate densities in each gas phase.
Table 3 lists column and number densities and Table 4 lists pillar masses.

####  5.1.1 13CO Column Density

We estimate the molecular gas column density assuming optically thin
13CO(J=1-0) emission following Tiwari et al. (2021) (their Appendix E) and
Mangum & Shirley (2015). We assume that the J=1-0 lines of 12CO and 13CO share
the same excitation temperature along each line of sight, that 12CO(J=1-0) is
optically thick everywhere, and that the beam filling factor is unity. We
adopt an isotopic ratio ${{}^{12}{\rm CO}}/{{}^{13}{\rm CO}}=^{12}{\rm
C}/^{13}{\rm C}=44.65$ based on the galactocentric radius-dependent expression
given by Yan et al. (2019) using a galactocentric radius of 6.64 kpc
calculated using Equation 2 by Brand & Blitz (1993). We convert 12CO column
density to molecular hydrogen column density $N({\rm H}_{2})$ using the
abundance ratio ${{}^{12}{\rm CO}}/{\rm H}_{2}=8.5\times 10^{-5}$ (Tielens,
2021). Hydrogen is assumed to be entirely in the molecular phase where CO is
detected, so that total hydrogen column density (defined as ${N_{\rm
H}}={N({\rm H})}+2{N({{\rm H_{2}}})}$) is ${N_{\rm H}}=2{N({{\rm H_{2}}})}$.
Gas mass is calculated by summing $N_{\rm H}$ over the boxes shown in Figure
12 and converting to mass using a mean molecular weight $\mu=1.33$.

We estimate line-of-sight-averaged densities by assuming cylindrical symmetry,
so that line-of-sight distance through the feature is equal to the angular
width of the features, and dividing column density by that distance. $H_{2}$
column densities and gas masses are listed in Tables 3 and 4, respectively.
Uncertainties are estimated by propagating the channel RMS noise of the 12CO
and 13CO(J=1-0) observations and, for conversion of column density into mass,
the heliocentric distance uncertainty, into the procedure described above.

Figure 12: The top-left, top-right, and bottom-left panels show hydrogen
column densities $N_{\rm H}$ derived from C+, CO(J=1-0), and FIR dust
emission. The bottom-right panel shows the JWST F335M image in greyscale as a
reference. All four plots are on the same size scale and the grey gridlines
are in the same place. On the bottom-right panel, blue boxes show the
integration areas for the pillar masses in Table 4, and red boxes show the
same for Table 6. Orange squares mark the locations from which column
densities are sampled in Table 3. In Table 3, from P1a-edge to Shared-Base-
Mid, locations are listed in order of decreasing declination along P1. The P2
locations are also in order of decreasing declination along P2.

####  5.1.2 C+ Column Density

We do not detect the [13C II] line towards the Pillars, so we estimate an
upper limit on the C+ column density using the [C II] channel noise $T_{\rm
rms}$ as a detection limit. We estimate an upper limit on the [13C II]
hyperfine component at +11.2 ${\rm km~{}s^{-1}}$ with respect to the [12C II]
line (see Table 1 in Guevara et al. 2020) since it is sufficiently far from
the [12C II] line center given the observed linewidths, and we assume that the
isotopic ratio ${{\rm{}^{12}C^{+}}}/{{\rm{}^{13}C^{+}}}=\ ^{12}{\rm
C}/^{13}{\rm C}$ and use the same value as in Section 5.1.1. Equation 4 in the
paper by Guevara et al. (2020) relates the observed line brightness of [12C
II] to [13C II] given the [12C II] optical depth. We find an upper limit of
$\tau_{12}\lesssim 1.3$ using the brightest [12C II] spectrum which lies
towards P1b. Bonne et al. (2023b) find a similar upper limit on [C II] optical
depth towards the DR21 ridge using the same method.

We use the method of estimating 12C+ column density without a detected [13C
II] line described by Okada et al. (2015), and we assume a constant ${{T_{\rm
ex}}}$ rather than a constant optical depth. The adopted ${{T_{\rm ex}}}=107$
K is calculated using the highest [C II] line brightness towards the Pillars
and our upper limit of $\tau\leq 1.3$, and ${{T_{\rm ex}}}$ is therefore a
lower limit towards that location (considering the emitted radiation
temperature $T_{\rm R}$ constant, $\tau$ and ${{T_{\rm ex}}}$ have an inverse
relationship in the optically thin and $\tau\sim 1$ regimes). Assuming a fixed
${{T_{\rm ex}}}$ across the field is equivalent to assuming that density and
kinetic temperature are fixed across the map, so that only the column density
influences the observed line brightness. Kinetic temperature is determined by
the density $n$ and radiation field $G_{0}$ according to PDR models. Both $n$
and $G_{0}$ may be slightly higher at the heads of the three Pillars which
contain dense molecular gas and are closer to the stars, but the column
densities there are not high enough to probe ${{T_{\rm ex}}}$. ${{T_{\rm
ex}}}$ must therefore be sampled from the higher column density P1b location.
We will overestimate column density towards dense, highly irradiated locations
like the pillar heads, and we will underestimate column density towards less
dense, poorly illuminated locations along pillar bodies. If ${{T_{\rm
ex}}}\sim 150$ K, column densities would be approximately halved around the
map, and it is unlikely that ${{T_{\rm ex}}}$ is considerably lower than 100
K. We convert $N({{\rm C^{+}}})$ to atomic hydrogen column density $N({\rm
H})$ using the abundance ratio ${\rm C}/{\rm H}=1.6\times 10^{-4}$ (Sofia et
al., 2004) and calculate column densities assuming a filling factor $\eta=1$.
We assume hydrogen to be entirely in the atomic phase where C+ is present so
that ${N_{\rm H}}={N({\rm H})}$, but we note that some [C II] emission may
originate from CO-dark molecular gas (Pabst et al., 2017).

The spatially variable [C II] background discussed in Appendix B contributes
to the total column density derived in the region. To determine just the
contribution to the Pillars requires subtracting a background. We calculate
column densities with the total intensities (no background spectrum
subtraction) and subtract out a column density background sampled from the
same locations as the spectral background. Within these background sample
regions, we take the mean column density to be the background and the standard
deviation to be a systematic uncertainty on both target and background column
density due to the variability of the background. Statistical uncertainty from
the RMS channel noise is propagated through the process described above,
though this source of uncertainty is reduced by the sum over pixels to obtain
mass. Background-corrected column densities are listed in Table 3. Estimates
of the systematic uncertainties and the 1-$\sigma$ statistical uncertainties
are given in the table’s caption.

We find the highest C+ column density towards P1b rather than P1a despite the
higher molecular gas column density towards P1a. Mass is calculated as
described in Section 5.1.1. The atomic hydrogen mass is typically $\sim$30% of
the molecular gas mass.

Table 3: Column and number densities.

|  |  | C+ |  | CO |  |  | Dust |
---|---|---|---|---|---|---|---|---|---
Location | Width |  |  $N_{\rm H}$ | $n_{{\rm H}}$ |  $N_{\rm H}$ | $\sigma_{\rm stat}$ | $n_{{{\rm H_{2}}}}$ |  $N_{\rm H}$ | $\sigma_{\rm stat}$
Name | (″) | (pc) | ($10^{21}$ cm-2) | (cm-3) | ($10^{21}$ cm-2) | ($10^{21}$ cm-2) | (cm-3) | ($10^{21}$ cm-2) | ($10^{21}$ cm-2)
P1a-edge | 17 | 0.14 | 20.9 | $4.7\times 10^{4}$ | 31 | 8 | $3.5\times 10^{4}$ | 59 | 12
P1a-center | 45 | 0.38 | 13.9 | $1.2\times 10^{4}$ | 286 | 21 | $1.2\times 10^{5}$ | 139 | 24
E-Thread | 20 | 0.17 | 8.2 | $1.6\times 10^{4}$ | 31 | 9 | $3.0\times 10^{4}$ | 22 | 7
W-Thread | 10 | 0.08 | 3.8 | $1.4\times 10^{4}$ | 29 | 9 | $5.5\times 10^{4}$ | 24 | 9
E-Horn | 13 | 0.11 | 14.8 | $4.4\times 10^{4}$ | 80 | 14 | $1.2\times 10^{5}$ | 22 | 6
W-Horn | 14 | 0.12 | 1.7 | $4.6\times 10^{3}$ | 125 | 15 | $1.7\times 10^{5}$ | 54 | 15
Shared-Base-E | 33 | 0.28 | 29.2 | $3.4\times 10^{4}$ | 59 | 13 | $3.4\times 10^{4}$ | 59 | 10
Shared-Base-Mid | 22 | 0.19 | 18.0 | $3.1\times 10^{4}$ |  $\leq$6 | 6 | … | 24 | 7
P2-head | 18 | 0.15 | 5.0 | $1.1\times 10^{4}$ | 114 | 14 | $1.2\times 10^{5}$ | 95 | 22
P2-clump | 9 | 0.08 | 6.3 | $2.7\times 10^{4}$ | 65 | 11 | $1.4\times 10^{5}$ | 59 | 16
P3-head | 13 | 0.11 | 3.6 | $1.1\times 10^{4}$ | 67 | 11 | $9.8\times 10^{4}$ | 58 | 22
Ridge | 13 | 0.11 | 19.2 | $5.7\times 10^{4}$ | … | … | … | 63 | 13

Note. — Column and number densities derived from C+, CO, and dust emission
measurements; the species is indicated above each column density column.
Column densities are sampled towards each location, shown in Figure 12, and
given in units of $10^{21}$ cm-2. The standard conversion between column
density and $A_{V}$ is $A_{V}={N_{\rm H}}/1.9\times 10^{21}~{}{\rm cm}^{-2}$
(Bohlin et al., 1978), so the column densities listed in the table can be
divided by 1.9 to estimate $A_{V}$. Column densities derived using the [C II]
line have statistical uncertainty $\sigma_{\rm stat}=0.6$ (table units) except
for location Shared-Base-E, which has $\sigma_{\rm stat}=0.7$, and an
additional systematic uncertainty of 1.5 or 1.7 (table units) for the northern
and southern background samples, respectively (see Appendices A and B), which
estimate the spatial variability of the C+ background (Section 5.1.2). Atomic
and molecular gas number densities $n_{{\rm H}}$ and $n_{{{\rm H_{2}}}}$ are
estimated using the C+ and CO measurements of column density by dividing
$N({\rm H})$ and $N({{\rm H_{2}}})$, respectively, by the assumed LOS width of
each feature. We assume ${N_{\rm H}}={N({\rm H})}$ in the atomic gas and
${N_{\rm H}}=2{N({{\rm H_{2}}})}$ in the molecular gas in order to estimate
density and pressure in each phase independently; however, some of the [C II]
emission could arise from CO-dark molecular gas. On-sky angular widths are
estimated using the JWST images and are adopted as the LOS widths through the
features assuming cylindrical symmetry.

Table 4: Atomic and molecular gas mass.

Feature | C+ |  |  | CO |  | dust |  |
---|---|---|---|---|---|---|---|---
Name | MH | $\sigma_{\rm stat}$ | $\sigma_{\rm sys}$ | M${}_{{\rm H}_{2}}$ | $\sigma_{\rm stat}$ | MH | $\sigma_{\rm tot}$ | Mtot
P1a | 24 | 1 | 7 | 79 | 3 | 83 | 23 | 103
P1b | 41 | 2 | 8 | 37 | 2 | 64 | 20 | 78
P2 | 34 | 2 | 7 | 70 | 3 | 95 | 33 | 103
P3 | 4.2 | 0.6 | 2.7 | 14 | 1 | 15 | 7 | 18

Note. — Gas mass $M$ derived for the pillars from [C II], 13CO(J=1-0), and
dust emission. Masses are given in $M_{\odot}$ and rounded; sums are
calculated from exact numbers. Statistical and systematic uncertainty
estimates are given where applicable. Total masses are considered to be the
sum of atomic gas mass from C+ and the molecular gas mass from CO.

####  5.1.3 Column Density from Dust Emission

We calculate column densities from FIR dust emission to compare to the
molecular and atomic gas column densities derived from CO and C+. Optical
depth at 160 µm is calculated directly using the 70 and 160 µm observations
and the method described by Tiwari et al. (2021). The $\geq$250 µm maps have
lower spatial resolution and do not resolve some Pillar features. We obtain
$\sim$14″ resolution in our final map by using only 70 and 160 µm. We subtract
the flux background in each band, sampled from the background regions
described in Appendix B (Figure 15), prior to calculation of dust properties
in order to isolate the Pillar emission. We divide the dust optical depth at
160 µm by the dust extinction cross section per H nucleus at 160 µm, $C_{\rm
ext,160}/{\rm H}$. We take the $R_{V}=3.1$ value of $C_{\rm ext,160}/{\rm
H}=1.9\times 10^{-25}$ cm-2/H from Draine (2003), following the method of
Tiwari et al. (2021). This yields $N_{\rm H}$, the column density of H nuclei
or total column density, directly. These column densities are, in theory,
equal to the sum of hydrogen column densities derived from C+, C, and CO as
dust is present at all phases, but they should be dominated by the relatively
larger molecular gas column. Column densities are listed in in Table 3 and
masses, calculated as described in Section 5.1.1, are listed in Table 4.

We measure dust temperatures ranging from 25 K in the cold clumps, such as the
pillar heads and the Western Horn, to 30 K along the lower column density
pillar bodies. Since we are using the 70 and 160 µm measurements for this
estimate, we are more sensitive to warm dust. Dust temperature varies along
the line of sight through the Pillars, with colder gas and dust embedded
within the dense pillar heads and other clumps, shielded from radiative
heating. Assuming a constant temperature along a line of sight with real
temperature variation can cause column density to be underestimated (Storm et
al., 2016; Howard et al., 2019). Including the 250 µm band increases
sensitivity to cold dust at the expense of spatial resolution and raises the
mass estimates by about 10% (these larger estimates are not included in the
table), which is within our uncertainty estimates described below.

The PACS handbook quotes a 5% calibration uncertainty on the 70 and 160 µm
measurements; this dominates the total uncertainty towards the pillar heads.
We estimate an uncertainty on our flux background subtraction in each band,
which is about 15% of the background and dominates total uncertainty outside
the bright pillar heads. The PACS observations include a statistical
uncertainty map; this is the smallest relative contribution to the total
uncertainty, but we include it anyway since including additional sources of
uncertainty is trivial using a Monte Carlo sampling method. Using such a
method, we estimate the total uncertainties on column density and mass by
computing $N_{\rm H}$ and $M_{\rm H}$ for 1000 realizations of the three
uncertainties added to the observations and adopting the 1-$\sigma$
uncertainties to be the difference between the 84th and 16th percentile values
divided by 2. To these uncertainties, we add in quadrature the 15% uncertainty
from the heliocentric distance (Kuhn et al., 2019). The median $N_{\rm H}$ and
$M_{\rm H}$ values are consistent with the values calculated without
considering uncertainty, so we record the latter values in the tables.

The dust-derived column densities tend to be closer to the CO-derived column
densities than to those from C+, satisfying our expectation that the dust
column is dominated by contribution from the molecular gas phase. The dust-
derived $N_{\rm H}$ values tend to be lower than the CO-derived values, with
the largest difference observed towards the bright center of P1a. This is
likely due to significant dust temperature variation along the line of sight,
which causes the average temperature to be higher and the optical depth to
then be underestimated. For this reason, we take the CO values as the better
gauge of molecular gas column density. Towards the northern edge of P1a, where
[C II] is brightest, the dust-derived $N_{\rm H}$ is larger than the CO-
derived value; we attribute this to the line of sight passing primarily
through warm atomic gas with less LOS temperature variation, as the sum of the
C+-derived and CO-derived column densities is consistent with the value from
dust.

###  5.2 Density

Line-of-sight averaged densities are derived by dividing the column density
measurements by estimated line-of-sight feature widths. We estimate the
angular size across each feature in Table 3 using the high-resolution JWST NIR
images and assume cylindrical symmetry to adopt the same width along the line
of sight. We use this width for both the atomic gas and the molecular gas,
even though they would be layered along lines of sight passing through
molecular gas. Densities $\sim 10^{5}$ cm-3 are derived from the CO- and dust-
based column densities. These densities are broadly consistent with the
observed HCO+, HCN, CS, and N2H+ brightness temperatures (listed in Appendix
A) according to estimates from the RADEX777http://var.sron.nl/radex/radex.php
radiative transfer software (van der Tak et al., 2007). Densities $\sim
10^{4}$ cm-3 are derived from C+-based column densities, which are consistent
with those estimated using the PDR Toolbox.

We adopt $n_{\rm H}={1.8\times 10^{4}}$ cm-3 as the atomic gas density by
taking the average $n_{\rm H}$ in Table 3. From this average we exclude the
Shared Base and Ridge positions, since their line-of-sight size depends
strongly on the adopted geometry which is very uncertain, and the Eastern
Horn, since we cannot isolate its [C II] emission from P1 emission behind it.
We adopt $n_{{{\rm H_{2}}}}={1.3\times 10^{5}}$ cm-3 as the molecular gas
density by taking the average $n_{{{\rm H_{2}}}}$ towards dense features
(heads, Horns, and the P2 clump) in Table 3. The molecular gas density we
derive is consistent with the estimates by Pound (1998), White et al. (1999),
and Levenson et al. (2000).

The critical density of [C II] in the presence of either H or H2 ranges from
3–${5\times 10^{3}}$ cm-3 at $\sim$100 K, so our adopted atomic gas density
implies that the kinetic temperature $T_{\rm K}$ should be similar to the
excitation temperature ${{T_{\rm ex}}}$. Recalling that our ${{T_{\rm ex}}}$
is a lower limit due to optical depth, we place a lower limit on atomic gas
$T_{\rm K}>100$ K. The derived molecular gas densities are greater than the
critical density of 12CO(J=1-0) by a factor of $\sim 50$, so the radiation
temperature $T_{\rm R}$, assumed equal to the observed $T_{\rm MB}$, of the
optically thick 12CO(J=1-0) line should be similar to the kinetic temperature
of the gas it probes. The kinetic temperature of the warm, outer layers of
molecular gas probed by 12CO must therefore be $T_{\rm K}\sim 50$–70 K for
most of the pillar bodies and reach $T_{\rm K}\sim 100$ K towards brightly
illuminated regions in P1a and P1b. The derived molecular gas density is close
to the critical densities $\sim 10^{5}$ cm-3 of HCO+ and HCN and slightly
below the critical densities $\sim{3\times 10^{5}}$ cm-3 of N2H+ and CS. The
HCO+ and HCN lines are optically thick according to RADEX given their assumed
abundance (Tielens, 2021) and the CO-derived column density $N({{\rm
H_{2}}})$, so their observed $T_{\rm MB}=T_{\rm R}\sim 15$–20 K likely traces
$T_{\rm K}\sim 20$ K molecular gas towards P1a.

###  5.3 Pressure Balance

We estimate pressures for the ionized, atomic, and molecular gas phases of gas
and comment on the possibility of pressure equilibrium. Pressures are listed
in Table 5. We follow White et al. (1999) in using the presence of pressure
equilibrium between the ionized and molecular gas as a proxy for whether or
not the ionization-driven shock has already been driven through the pillar
heads.

####  5.3.1 Ionized Gas

Our observations do not probe directly the physical conditions of the ionized
hydrogen (H II) phase, so we take values $n_{\rm e}\sim 1120$–$1800$ cm-3 and
$T_{\rm e}\approx 8000$–10,000 K from optical ion line studies of the pillar
surface by García-Rojas et al. (2006) and McLeod et al. (2015). Levenson et
al. (2000) note that temperature may be higher near the PDR-H II interface due
to photoelectric heating, a harder radiation field, and high density. We
estimate thermal pressure in the ionized gas flowing from the surface of the
pillar $P_{\rm i}/k_{\rm B}=2n_{\rm e}T=1.8$–$3.6\times 10^{7}$ K cm-3. We
neglect any sources of non-thermal pressure in the ionized gas near the pillar
surface; magnetic fields should not be important in the ionized phase due to
low densities and, while a turbulent linewidth was detected in the ambient H
II region (far from the pillar surface, where the photoevaporative flow does
not dominate gas kinematics) in M16 by Higgs et al. (1979), the radial
striations evident in the optical images presented by Hester et al. (1996)
indicate that the photoevaporative flow near the pillar surface is not
turbulent.

####  5.3.2 Molecular Gas

Pattle et al. (2018) discuss the pressures originating within the molecular
gas phase. We follow their discussion and recalculate pressures using our
observations of the molecular gas. We determined $n_{{\rm H_{2}}}={1.3\times
10^{5}}$ cm-3 in Section 5.2. Molecular gas temperature should be coupled to
dust temperature $T_{d}$ at densities $n\gtrsim 10^{5}$ cm-3 (White et al.,
1999). We find $T_{\rm d}\sim 25$ K towards the dense gas (Section 5.1.3),
comparable to the value $T\sim 20$ K found by White et al. (1999) using
350–2000 µm photometry, which is more sensitive to cold dust but less
sensitive to warm dust than the 70 and 160 µm photometry used for our
estimate. Thermal pressure in the molecular gas is $P_{{{\rm H_{2}}},\,{\rm
therm}}/k_{\rm B}=n_{{\rm H_{2}}}T\approx{3.2\times 10^{6}}$ K cm-3. We
estimate the pressure from turbulent support using the molecular line velocity
dispersion $\sigma_{\rm obs}\approx 0.6$ ${\rm km~{}s^{-1}}$ (thermal velocity
dispersion, $\sigma_{\rm therm}\sim 0.1~{}{\,{\rm km~{}s^{-1}}}$, contributes
less than 2% of $\sigma_{\rm obs}$) to be $P_{{{\rm H_{2}}},\,{\rm
turb}}/k_{\rm B}=\rho\sigma^{2}=(2n_{{\rm H_{2}}}\mu m_{\rm
H})\sigma^{2}\approx{1.5\times 10^{7}}$ K cm-3 ($\mu=1.33$), consistent with
pressure derived from the line widths observed by White et al. (1999).

Pattle et al. (2018) estimate the plane-of-sky magnetic field strength
$B=170$–320 $\mu$G within the Pillars, calculated using the Davis-
Chandrasekhar-Fermi method (Davis, 1951; Chandrasekhar & Fermi, 1953) and
therefore assuming that the turbulent velocity field is isotropic, which
implies magnetic pressure $P_{{{\rm H_{2}}},\,{\rm B}}/k_{\rm B}\approx
0.9$–$3.0\times 10^{7}$ K cm-3. The factor of 3 range in their pressure value
seems to be driven mostly by variation in the molecular line widths observed
by White et al. (1999); our observations indicate that some of this line width
can be attributed to line-of-sight confusion rather than turbulent velocity
dispersion and we observe single-component line widths near the lower end of
their range. We use the Equation 1 in the paper by Pattle et al. (2018) to
estimate $B_{\perp}\sim 320~{}\mu$G and $P_{{{\rm H_{2}}},\,{\rm B}}/k_{\rm
B}\approx{3.0\times 10^{7}}$ K cm-3 from our observed molecular gas velocity
dispersion $\sigma_{v}\approx 0.6$ ${\rm km~{}s^{-1}}$ along lines of sight
with limited line-of-sight confusion. The magnetically-dominated total
pressure within the molecular gas is $P_{{\rm H_{2}}}/k_{\rm
B}\approx{4.8\times 10^{7}}$ K cm-3.

The total pressure in the molecular gas is slightly higher than the pressure
in the ionized gas. The molecular and ionized gas phases are likely in
pressure equilibrium, with the extra pressure supporting the molecular gas
against self-gravity. We assume that the ionization-driven shock has already
passed through the molecular gas.

We observe a velocity gradient across the head of P2, which could indicate an
additional mode of rotational support in the pillar head. For comparison, we
estimate a rotational “pressure” $P_{{{\rm H_{2}}},\,{\rm
rot}}=\rho\omega^{2}r^{2}$ where $\omega r\approx 0.5$ ${\rm km~{}s^{-1}}$
based on the CO line velocities (see the PV diagram in Figure 9), finding
$P_{{{\rm H_{2}}},\,{\rm rot}}/k_{\rm B}\approx 10^{7}$ K cm-3 (not included
in our total pressure). This is comparable to turbulent pressure and less than
magnetic pressure (${3\times 10^{7}}$ K cm-3, see above), and could contribute
$\sim$20% of the total support to the pillar head, but we would expect the
head to be extended in the plane of rotation across the pillar if that were
so. Since we do not observe oblateness, we find it unlikely that rotational
support is a significant support mechanism in the pillar heads. Pillar
rotation is further discussed by Sofue (2020).

####  5.3.3 Atomic Gas

We determined the density of the warm, atomic gas to be $n_{{\rm
H}}\approx{1.8\times 10^{4}}$ cm-3. The lower limit on the ${{T_{\rm ex}}}$ of
the [C II] line indicates a kinetic temperature $T_{\rm K}>100$ K (Section
5.2), and temperatures up to $250$ K are consistent with the approximate $n$
and $G_{0}$ in the PDR models (Pound & Wolfire, 2023). This results in a
thermal atomic gas pressure $P_{{\rm H}}/k_{\rm B}=n_{{\rm H}}T\sim
2$–$5\times 10^{6}$ K cm-3. The atomic gas phase is likely hydrostatically
supported by turbulence similar to the molecular gas. We remove the thermal
velocity dispersion of [C II] $\sigma_{\rm therm}\sim 0.5~{}{\,{\rm
km~{}s^{-1}}}$ given the assumed $T_{\rm K}$ from the observed [C II] velocity
dispersions $\sigma_{\rm obs}\sim 1.1$–1.4 ${\rm km~{}s^{-1}}$ to find
$\sigma_{\rm turb}=(\sigma_{\rm obs}^{2}-\sigma_{\rm therm}^{2})^{1/2}\sim
1.0$–1.3 ${\rm km~{}s^{-1}}$. Adopting the turbulent velocity dispersion
$\sigma_{\rm turb}\approx 1.3$ yields a turbulent pressure of $P_{{\rm
H},\,{\rm turb}}/k_{\rm B}\sim{5\times 10^{6}}$ K cm-3 in the atomic gas. This
is an upper limit since the photoevaporative flow velocity $\sim 0.5$ ${\rm
km~{}s^{-1}}$ projected along the line of sight should contribute to the
observed width of the [C II] line, but we cannot disentangle these
contributions in our observations. We discuss the impact of this flow velocity
on the photoevaporative lifetime of the Pillars in Section 6.2. Ultimately
even the highest estimate of turbulent pressure in the atomic gas allowed by
the observations is an order of magnitude smaller than the total pressures in
the ionized and molecular phases in Table 5.

It is likely that ionization-driven shocks have already passed through the
Pillars, so the atomic gas should be in pressure equilibrium with the ionized
and molecular gas phases surrounding it. However, our estimate of the total
atomic gas pressure falls significantly short of the ionized gas pressure by
$P/k_{\rm B}\sim{1\times 10^{7}}$ K cm-3. The missing pressure is likely from
the magnetic field in the atomic gas whose strength must be $B\sim 200\mu$G.
Pellegrini et al. (2007, 2009) estimate comparable magnetic field strengths in
the PDRs in M17 and the Orion Bar and determine that the magnetic field is
important to the structure, particularly the width, of the PDR. We don’t have
an independent estimate of the field strength in the atomic layer, but scaling
the molecular phase $B$ field by the ratio of the molecular to atomic gas
densities gives a field strength of $B_{\perp}\sim 85~{}\mu$G and a
corresponding $P_{\rm B}/k_{\rm B}\sim{2\times 10^{6}}$ K cm-3, which is
insufficient for pressure equilibrium.

We conclude that ambipolar diffusion, the flow of neutrals past ions, has
weakened the magnetic field in the molecular gas faster than in the atomic gas
because the molecular gas has a lower ionization fraction. The pillar heads
are still magnetically supported now, but the magnetic field in the molecular
gas must have been larger in the past. The ambipolar diffusion timescale
$\tau_{\rm AD}\sim 10^{6}$ yr in the molecular gas is similar to the age of
the cluster, as we demonstrate in the next section. Carbon is ionized in the
atomic gas, so the ionization fraction there is $\chi_{\rm ion}\sim 10^{-4}$
and $\tau_{\rm AD}\sim 10^{9}$ yr.

###  5.4 Ambipolar Diffusion and Molecular Cloud Collapse

Mouschovias (1991) determined that the reduction in magnetic support due to
ambipolar diffusion could precipitate the collapse of molecular cores into
stars. We estimate the possibility of such a collapse in the Pillars before
they are photoevaporated. Critical mass under magnetic field support can be
expressed as

$M_{\rm c,B}=9.7\Bigg{(}\frac{R}{0.1~{}{\rm
pc}}\Bigg{)}^{2}\Bigg{(}\frac{B}{100~{}\mu{\rm G}}\Bigg{)}M_{\odot}$ (2)

The magnetic critical masses, calculated using the assumed radii in Table 6,
for the heads of P1a, P2 and P3 are 83, 21, and 10 $M_{\odot}$, respectively,
while the estimated H2 masses are 72, 14, and 6 $M_{\odot}$. For all the
pillar heads, but particularly P1a whose geometry is complex (see Section
3.2), the $R^{2}$ dependence carries significant uncertainty into the critical
mass determination from our assumptions about line-of-sight width, but we can
conclude that all the heads are nearly critical. Lefloch & Lazareff (1994)
present calculations to determine the gravitational stability of magnetically
supported globules in Section 5.9 of their paper, and these suggest that the
pillar heads are unstable given their current magnetic field strengths.

For magnetically subcritical clouds, gravitational collapse is controlled by
the process of ambipolar diffusion of the supporting magnetic field. The
ambipolar diffusion timescale is

$\tau_{\rm AD}\approx{2\times 10^{5}}\Big{(}\frac{\chi_{\rm
ion}}{10^{-8}}\Big{)}~{}{\rm yr}$ (3) $\chi_{\rm ion}\approx{2\times
10^{-7}}\Bigg{(}\frac{10^{4}~{}{\rm cm}^{-3}}{n_{{\rm H_{2}}}}\Bigg{)}^{1/2}$
(4)

where the degree of ionization $\chi_{\rm ion}$ is calculated for a typical
cosmic ray ionization rate $\zeta_{\rm CR}={3\times 10^{-17}}$ s-1 (Elmegreen,
1979). For the adopted $n_{{\rm H_{2}}}={1.3\times 10^{5}}$ cm-3, we find
$\chi_{\rm ion}\sim{6\times 10^{-8}}$ and $\tau_{\rm AD}\sim 1$ Myr which is
similar to both the age of the cluster and Pillars as well as the estimated
photoevaporative lifetime of the Pillars (Section 6.2). The models by Bergin
et al. (1999), which use physical conditions similar to those observed towards
the pillar heads, indicate that $\chi_{\rm ion}<10^{-8}$ is more typical of
these massive cores; this would mean ambipolar diffusion weakens the magnetic
field more quickly than expected. Adopting this lower ionization fraction
would yield $\tau_{\rm AD}\sim 0.2$ Myr, smaller than both the age of the
system and the photoevaporative lifetime. It is possible that any of the
pillar heads may collapse to form stars before they are photoionized as
ambipolar diffusion weakens the magnetic fields supporting them against
gravitational collapse. Protostars have been observed towards the heads of P1a
and P2, so fragmentation and local collapse may have already taken place
(Sugitani et al., 2002; Indebetouw et al., 2007; Linsky et al., 2007).

The interplay between gravity, magnetic fields, and turbulence in context of
cloud collapse is a complicated one which we simplify in our discussion here;
see the review by Hennebelle & Inutsuka (2019). The ambipolar diffusion
timescale is slow compared to the dynamic timescales of clouds. Quicker
alternatives for overcoming magnetic support are available, such as turbulence
or other methods of locally enhancing the mass-to-magnetic-flux ratio
(Vázquez-Semadeni et al., 2011; Bailey & Basu, 2014). If magnetic support were
overcome more quickly, the magnetically supported and critical Pillar heads
would be even more likely to fragment and collapse into stars before they
evaporate.

Table 5: Pressures in the molecular, atomic, and ionized gas phases.

Component | $P_{therm}$ | $P_{turb}$ | $P_{B}$ | $P_{tot}$
---|---|---|---|---
Molecular gas | 3 | 15 | 30 | 48
Atomic PDR gas | 2–5 | 5 | 2* | 9–12
Ionized flow | 18–36 |  |  | 18–36

Note. — Pressures in the molecular, atomic, and ionized gas near P1a.
Pressures are expressed as $P/k_{\rm B}$ in units of ($10^{6}$ K cm-3). Total
pressures are calculated by summing across lower and upper estimates,
respectively. The value ranges are intended as approximate estimates of our
uncertainty in the pressures due to uncertain, or multiple, estimates of
density or other properties. The atomic gas magnetic pressure (*) in this
table is calculated by scaling down the magnetic field in the molecular gas by
the ratio of mass densities in the two phases, but we suspect it is much
larger as discussed in the text.

###  5.5 Constraints on Pillar Age

Simulations by Williams et al. (2001) indicate that the shock can pass through
the dense, molecular pillar head in $10^{5}$ yr while the ionization front
takes $\sim 0.5\times 10^{6}$ yr to photoevaporate the pillar, comparable to
the $\sim 10^{6}$ yr age of the system. Our observations are consistent with a
shocked, equilibriated pressure structure, in which the shock passed through
long ago and left the ionization front behind. This relieves us of the tight
upper limit of $\sim 10^{5}$ yr on the pillar age required for unshocked
pillar heads (White et al., 1999; Williams et al., 2001). Our estimated
photoevaporation timescales $\sim 10^{6}$ yr (Section 6.2) and the NGC 6611
cluster age $\sim{2\times 10^{6}}$ yr also favor the post-shock, $\sim$106-yr-
old Pillars scenario. See also the discussion of timescale tension by Williams
et al. (2001).

##  6 Discussion

###  6.1 Warm Gas Dynamics

Along lines of sight towards the Shared Base and the head of P2, the peak [C
II], [O I], and CO line velocities are blueshifted by $\sim$0.5 ${\rm
km~{}s^{-1}}$ with respect to the denser molecular gas tracers. The [C II]
line, but not the CO, is similarly blueshifted towards the head of P3. This is
inconsistent with a photoevaporative flow from an illuminated surface that is
mostly facing away from us, but is consistent with gas flowing down these
pillars away from the star (towards us in projection along the LOS).

Simulations by Lefloch & Lazareff (1994) predict that less dense, outer layers
of gas move down the pillar bodies more quickly than denser gas (see their
Figure 4c–f). These simulations also predict blue line wings towards the head
of a globule (oriented away from us, as P1b and P2 are) at two different
evolutionary stages: a wing due to compressed gas on the far side of the head
moving towards the observer during the first $\sim$0.1 Myr while the cloud is
collapsing into a globule, and a wing during the stable cometary phase at an
age of $\sim$1 Myr as small clumps are ejected and accelerated along the
globule. As we mention in Section 3.5, the wishbone tails of P3 resemble the
“ears” which emerge in these simulations; comparisons could be drawn between
the later stages of these ears and the Threads. Observational comparisons
could be made between the Threads, ears, and tails and the “legs” observed by
Bonne et al. (2023a) towards IC 63.

Since [C II], [O I], and 12CO should all trace warmer, outer layers of gas,
the observed velocities can be explained if these layers are moving faster
down the pillar than the colder, denser gas deeper in the pillar heads. We, in
our [C II] and molecular lines, and Pound (1998) in CO and McLeod et al.
(2015) in optical ion lines, have observed radial velocity gradients along the
Pillars which are thought to trace gas accelerated away from the star as it
flows down the pillar. Flow of gas away from the head down the body is
predicted in simulations of certain pillar configurations by Mackey & Lim
(2010) and of cometary globules by Lefloch & Lazareff (1994). To explain our
blueshifted outer layers gas, we need only require that the phenomenon acts
more strongly on less dense gas, accelerating it more quickly near the head so
that we see a relative velocity shift there but not further down the pillar
body where all phases have been accelerated. The driving force is unknown, but
we suggest a few possibilities. This could be a shearing flow generated by the
action of stellar photons on the slanted pillar structure. The action is
transmitted in the warm PDR surface layers, so they respond first. Deeper
self-gravitating layers will be slowly carried along through shearing action.
Other candidates for the driving force could be the rocket effect or shocks
driven by the ionization front, perhaps in combination with surface geometry.
Stellar winds, while a good candidate for transferring momentum radially away
from the exciting stars, do not reach the surface of the pillars. The
striations detected in optical images by Hester et al. (1996) indicate that
photoevaporated material flows away from the ionization front uninterrupted
for at least a few tenths of a parsec, so the winds must terminate against the
H II region further away towards the stars (see the diagram in the left panel
of Figure 19 in the paper by Westmoquette et al. 2013).

Towards the Eastern Thread, the [C II] and 12CO lines are blueshifted with
respect to CS and other molecules. P1a faces towards, not away from, the
observer, so this velocity shift is not well explained by these mechanisms.
The velocity shift towards the Eastern Thread may be due to an entirely
different phenomenon than the velocity shifts towards the Shared Base, P2, and
P3.

###  6.2 Pillar Lifetimes

We estimate the photoevaporative lifetime of the pillars based on our
observations. We adopt the equation for mass loss rate
$\dot{M}=Av_{f}\rho_{\rm H}$ from Gorti & Hollenbach (2002), where $A$ is the
total evaporative surface area, $v_{f}$ is the photoevaporative flow velocity
from the surface, and $\rho_{\rm H}$ is the mass density at the base of the
flow.

Table 6: Pillar head and neck masses, mass loss rates, and photoevaporative
lifetimes.

Feature | R(, H) | R(, H) | A | $M_{H}$ | $M_{H_{2}}$ | $M_{\rm tot}$ | $\dot{M}$ | $t$
---|---|---|---|---|---|---|---|---
Name | (″) | (pc) | (pc-2) | ($M_{\odot}$) | ($M_{\odot}$) | ($M_{\odot}$) | ($M_{\odot}$ Myr-1) | (Myr)
P1a Head | 20 | 0.17 | 0.18 | 15 | 72 | 86 | 54 | 1.6
P1a Neck | 20, 38 | 0.17, 0.32 | 0.34 | 8 | 7 | 15 | 103 | 0.1
P2 Head | 10 | 0.08 | 0.04 | 2 | 14 | 16 | 14 | 1.2
P2 Neck | 10, 28 | 0.08, 0.24 | 0.13 | 3 | 21 | 24 | 38 | 0.6
P3 Head | 7 | 0.06 | 0.02 | 1 | 6 | 7 | 7 | 1.0

Note. — Pillar heads are modeled as hemispheres and necks as cylinders; the
R(, H) columns give radius and, if applicable, height. The adopted height for
each pillar neck matches the height of the box used to sum over column
density. Mass values are rounded to the nearest whole number. Adopted atomic
gas density at the base of the photoevaporative flow is $n_{\rm H}={1.8\times
10^{4}}$ cm-3. Adopted flow velocity is 0.5 ${\rm km~{}s^{-1}}$.

Warm atomic gas should flow along the pressure gradient through the PDR
towards the ionization front at a $v_{f}$ comparable to, but no greater than,
the sound speed in the atomic gas $c_{\rm H}=(P_{{\rm H},therm}/\rho_{\rm
H})^{1/2}\sim 1~{}{\,{\rm km~{}s^{-1}}}$. The observed [C II] velocity
dispersion gives an upper limit $v_{f}\lesssim\sigma_{\rm obs}\approx
1.4~{}{\,{\rm km~{}s^{-1}}}$, though we attribute most of the observed
velocity dispersion to turbulence (Section 5.3). We adopt $v_{f}\sim
0.5~{}{\,{\rm km~{}s^{-1}}}$.

Mass density at the base of the flow $\rho_{\rm H}=\mu m_{\rm H}n_{\rm H}$,
where we use $\mu=1.33$ and the atomic gas density $n_{\rm H}={1.8\times
10^{4}}$ cm-3 from Section 5.2. Evaporative surface area $A$ is estimated
based on the on-sky sizes of the Pillars in high-resolution optical and NIR
images. In order to better understand the evolution of the Pillars, we
estimate the evaporative lifetimes of the heads and bodies separately.
Evaporation of low column density regions towards the pillar bodies is of
particular interest, since we detect significant [C II] line and PAH emission
from these regions. We define the pillar “necks” as the lower column density
segments of the body between the head and any dense clumps along the body. We
use the lack of detected N2H+ emission to distinguish the neck from the head
and rest of the body. For P1a, the neck comprises both Threads, and for P2 the
neck extends between the pillar head and the clump. We do not resolve P3 well
enough to define a neck. The boxes used to integrate the head and neck masses
are shown in the rightmost panel of Figure 12. We model pillar heads as
hemispheres and necks as cylinders, and we select a single typical radius for
both components of each pillar based on the angular sizes in Table 3. The
height of the neck is the length of the box used to integrate its mass. The
adopted dimensions and corresponding mass loss rates are listed in Table 6.
The surface area is the only variable changed from region to region, as the
number densities in Table 3 do not vary by more than a factor of a few and the
total [C II] velocity dispersion only varies by $\sim$50%. Column densities
calculated using C+ and CO are summed pixel-by-pixel within the boxes to find
the atomic and molecular gas masses, and the total mass of each feature is
considered to be the sum of the mass in both phases. The photoevaporative
lifetime of each feature is the mass divided by the mass loss rate. The masses
and photoevaporative lifetimes are listed in Table 6.

The pillar heads will evaporate in $\sim$1–2 Myr by our estimate, consistent
with the estimate by McLeod et al. (2015) of $\sim$3 Myr for these pillars as
well as the estimate by Westmoquette et al. (2013) of $\sim$2 Myr for the
pillars surrounding NGC 3603. This is similar to the $\sim$2 Myr age of the
NGC 6611 system, so the Pillars likely formed with roughly a few times as much
mass as they have now and have photoevaporated over the last $\sim$Myr down to
their present-day masses. As we discussed in Section 5.3, the photoevaporative
lifetime is at least as large as $\tau_{AD}$, so the heads of the Pillars may
still collapse to form stars before they are completely photoevaporated.

The pillar necks appear to evaporate on a shorter timescale of $\sim$0.1 Myr
due to their much larger surface areas. If we have overestimated their mass
loss rates, the culprit could be a diminished flow velocity which is obscured
in the observed velocity dispersion by the more dominant turbulent velocity
dispersion. In the case of Pillar 1, the flow velocity would have to be
diminished by a factor of $\sim$10\. We do observe thinner (by $\lesssim$0.5
${\rm km~{}s^{-1}}$) line profiles towards the Threads in P1a and the neck of
P2, so this order of magnitude decrement in flow velocity between the head and
the neck is possible and must be observationally confirmed with, for example,
velocity resolved radio recombination lines or optical ion lines. $G_{0}$,
measured from FUV emission, is diminished only by $\sim$60% from the Head to
the Neck of P1a, and by $\sim$30% in P2.

We briefly consider the possibility that the necks do evaporate $\sim$10x more
quickly: either they will evaporate and cause the dense clumps towards the
heads to appear as cometary globules like those observed towards Cygnus X
(Schneider et al., 2016) or at the 0.01 pc scale elsewhere in M16, or they are
being replenished by bulk gas flows from the head like those discussed in
Section 6.1, which may hasten the dissipation of the pillar heads by a factor
of $\sim$2\. The meaningful difference between the two evolutionary scenarios
is whether the Pillars are observed as globules or coherent pillars for the
last $\sim$1 Myr of their lives; the work of Schneider et al. (2016) suggests
the former.

##  7 Conclusion

We observed the Pillars of Creation in M16 in velocity-resolved [C II] and [O
I] for the first time and detected PDRs which are spatially and kinematically
associated with the molecular gas. We combine these data with observations of
APEX, CARMA, BIMA, Herschel, Spitzer, and JWST to obtain a multi-layer view of
these gas structures.