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

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 μm line. These observations are velocity resolved to 0.5 km s−1 and are analyzed alongside a collection of complimentary data with similar spatial and spectral resolution: the [O i] 63 μ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 nH2∼1.3×105 cm−3 and an average atomic gas density n H ∼ 1.8 × 104 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 ⊙ for P1a, P1b, P2, and P3, respectively, and evaporation times of ∼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.

The H II region associated with M 16 is powered by the early-type stars in the NGC 6611 cluster, a ∼ 2×10 4 M ⊙ cluster whose most massive member is estimated to be ∼ 80 M ⊙ (Hillenbrand et al. 1993;Hill et al. 2012).The cluster, at a heliocentric distance of ∼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. (2007Guarcello et al. ( , 2009Guarcello et al. ( , 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 ∼1 Myr-old stars, broadly consistent with the ∼2 Myr age determination of NGC 6611 by other authors, whose formation event may have been triggered by a supernova shell's arrival ∼3 Myr ago.We adopt the distance of 1740 ± 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(Belikov et al. , 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 ∼ 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, ∼ 10 5 yr, may be shorter than their estimated photoevaporation timescale ∼ 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 welldefined.We use the JWST images throughout this study to measure the on-sky angular widths of features.

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 andWolfire 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 ∼ 100 K) than the molecular gas (T ∼ 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 ≳ 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. (2012Schneider et al. ( , 2021) ) towards pillars and globules in Cygnus X.
At A V ≳ 5, carbon is found primarily in the molecular phase as carbon monoxide (CO).From A V ∼ 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 12 CO(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 ≳ 10.
Molecular lines more sensitive to dense gas (n ≳ 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 velocityresolved [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.  (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 ′ ≈ 0.5 pc at 1740 pc line of sight.
M16 was observed with upGREAT 1 on SOFIA on 9 flights between 2019 and 2022 from Palmdale, California and Tahiti in the [C II] 3 P 3/2 → 3 P 1/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 ∼ 590 arcmin2 was mapped in the on-thefly (OTF) mode and atmospheric calibration was done with the GREAT pipeline (Guan et al. 2012).A Fast Fourier Transform Spectrometer (FFTS) with 4 GHz 1 upGREAT 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.
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 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 MB using an average main beam efficiency of 0.65.The forward efficiency is η 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) 2 .
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 µ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 km s −1 .

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 H 2 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 Darray configuration track was obtained on November 12, 2012.The antenna signals were transmitted to the eightband spectral line correlator.Four bands were used to observe the spectral lines HCN(J=1-0), HCO + (J=1-0), N 2 H + (J=1-0), and CS(J=2-1) in 7.8 MHz bandwidths with spectral resolution 24 kHz/channel (∆V ∼0.08 km s −1 ).Four bands were used to measure continuum at ∼ 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 ∼ 10 mJy; it is not discussed in this study.)The Eagle was observed with a 37-point hexagonal mosaic centered on (α, δ) J2000 = (18 h 18 m 51.s 29, −13 • 15 ′ 02.32 ′′ ).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), 13 CO(J=1-0) and C 18 O(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(Sault et al. , 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 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 N 2 H + (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 ∼1 km s −1 from each other.Since these are not well separated given the ∼1 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 N 2 H + (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.
The array is arranged in a hexagonal configuration around a central pixel with a spacing of about two beam widths (θ MB = 18.2 ′′ 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 × 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 (α, δ) J2000 = (18 h 20 m 46.s 3, −13 • 14 ′ 56 ′′ ).
The mapped region, centered at 18 h 18 m 35.s 7, −13 • 43 ′ 31.0 ′′ and of size 30 ′ × 22 ′ at −40 deg angle (CCW against positive RA), was split into 2 × 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 T MB (main-beam efficiency η MB = 0.68 at 345.8 GHz).All linear baselines were removed and all data re-sampled to 0.1 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.9mm) in the 12 CO(J=6-5) transition using the SEPIA660 receiver and the FFTS1 backend spectrometer on the APEX telescope.The instrument was tuned so that the 12 CO(J=6-5) line rest frequency 691.473GHz lay in the upper sideband, whose bandwidth is 4 GHz and frequency resolution is 61 kHz.The OTF-mapped region is a 4.3 ′ × 4.3 ′ square centered at 18 h 18 m 51.s 2, −13 • 50 ′ 0 ′′ .Spectra are calibrated in T MB (main-beam efficiency η 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 km s −1 .

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 website 4 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).
4 https://irsa.ipac.caltech.edu/data/SPITZER/GLIMPSE/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.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-α 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 ∼ 200 AU resolution.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 α, δ = 18 h 18 m 56.s , −13 • 51 ′ 20 ′′ in Figure 3), while the CS integrated intensity map overlaid in that figure unveils the density structure beneath.
Pillar 1 appears between V LSR = 24-26 km s −1 in the channel maps in Figure 4 and Pillars 2 and 3 appear between V LSR = 20-23 km s −1 .The ∼10 pc scale filaments observed by Hill et al. (2012) and Xu et al. (2019) lie at similar velocities V LSR = 20, 22.5, and 25 km s −1 .Much of the [C II] emission in the channel maps is diffuse emission surrounding the Pillars between V LSR ∼ 24-27 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 V LSR ∼ 25 km s −1 [C II] component.The pillar system connects towards the southeast to a ∼10 pc scale sys- 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 rootmean-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.
tem of emission features (only a small part of which is shown in Figure 4) at V LSR ∼ 25 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 V LSR ≥ 29 km s −1 reveal a faint, ring-like feature ∼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.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.

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.] are integrated between VLSR = 18-27 km s −1 , as they tend to have longer blue wings, and molecular lines not otherwise specified are integrated between VLSR = 19-27 km s −1 .N2H + and HCN are integrated between VLSR = 12.6-32 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 km s −1 for C 18 O(J=1-0); [4, 24, 44] K km s −1 for 13 CO(J=3-2); and [9, 18, 27, 45] ×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.
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  [25,50,75,100,125,150,175,200,225] K 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 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.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.01pc 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.

P1a: Threads and Cap
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 ∼0.3 pc.There are red-to-blue radial velocity gradients of order −1 km s −1 pc −1 along each Thread towards the pillar head (see between 70-100 ′′ and V LSR < 27 km s −1 in the position-velocity diagram along P1 in Figure 8).The Eastern Thread is spatially broader, ∼0.15 pc wide in CO emission just below the head of the pillar, and is redshifted by ∼1 km s −1 with respect to the thinner (∼0.07 pc) Western Thread.The Eastern Thread has a transverse gradient in radial velocity, evident in the 12 CO 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 N 2 H + line emission towards P1a around V LSR = 23-26 km s −1 extends slightly towards the northernmost part of the Western Thread, indicating higher density.
The Cap is a lower velocity (V LSR = 22-24 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  A. Wolpa (NOIR-Lab/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.net a (Lang et al. 2010).The colors show the O III line at 499 nm (blue), the Hα line at 656 nm (green), and the S II line at 672 nm (red).a https://astrometry.netP1a'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 N 2 H + 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 12 CO lines, and is present but significantly diminished in HCO + when compared to CS (see the spectra in Appendix A).

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 ∼10 ′′ .These ∼0.1 pc diameter clumps are particularly pronounced in the molecular line channel maps, indicating the presence of dense molecular gas, and the N 2 H + and C 18 O lines are faintly detected towards the Western Horn but not the Eastern Horn.The Western Horn is brighter in longer wavelength (> 160 µ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 pillarlike 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 (∼ −3 km s −1 pc −1 ) smoothly connecting P1 to P2.The [C II] and 12 CO lines contain excess low-velocity emission between V LSR ≈ 21-22 km s −1 compared to higher critical density molecular lines like HCO + ; we show example spectra in Appendix A.

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, ∼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 V LSR = 22-23 km s −1 towards the pillar head; high density tracers like N 2 H + peak closer to V LSR = 23 km s −1 while warmer, lower density tracers like [C II] and 12 CO(J=1-0) are blueshifted by about 1 km s −1 and all 12 CO 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 C 18 O and N 2 H + 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 13 CO 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.

P3
Pillar 3, the smallest of the three main Pillars, presents in our molecular line observations as a ∼0.38 pc long pillar oriented roughly parallel to P1 and P2 with two tails which extend ∼0.34 pc in either direction at ∼100 • angles from the body in a "wishbone" shape also seen in optical images.The head of the pillar is brightest at V LSR ≈ 21.2 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 ≲ 0.3 km s −1 through the head.C 18 O and N 2 H + 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.

Ridge
The Ridge and P4 (Section 3.7) are detected in most molecular lines and are each associated with a compact clump of N 2 H + 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 N 2 H + half-power beam.
The Ridge, lying around 0.5 pc southeast of the main Pillars' tails and appearing in the same V LSR ≈ 25 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 ∼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 V LSR = 24-27 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 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, ∼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 12 CO(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.

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 (V LSR ∼ 21 km s −1 ) up its western side to its point (∼ 23 km s −1 ) and back down its northeastern corner (∼ 25 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.

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-α 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 ∼0.5 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 • , 40 • , and 40 • 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 ≳0.5 pc (see Sections 4.1.2and 5.1), and our analysis of the FUV radiation field places the heads of P1a and P2 ∼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-α, and the Ridge farther from the observer based on its relative brightness in Pa-α (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.

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 ∼1 km s −1 rather than from ionized gas which would have a turbulence-dominated line width closer to 10 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 ∼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.

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, 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-ofsight 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 ∼1 pc along the line of sight; see Section 3.8 for more detail.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.N 2 H + 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.

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 V LSR = 23-25 km s −1 .The bases of the Horns are connected by warm, low column density 12 CO-bright gas between V LSR = 23-24 km s −1 (White et al. 1999).The Shared Base lies at V LSR = 22-23 km s −1 .The Shared Base and the Horns are spatially and dynamically connected but are separated by about 0.13 pc and 1 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 (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 ∼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 ≳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 ≥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.

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 V LSR = 20-22 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.

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 (∼14 ′′ , 0.1 pc) is comparable to their beams.At 3.3 µm, the head of P3 is capped by a thin (≈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.

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 (≈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 ≳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.

P4
The two sides of the angular P4 both host edge-on PDRs at their surfaces.We assume the low 1.87 µm Pa-α 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.

Sources of Illumination
The brightest members of NGC 6611 are ∼ 2.5 pc from the Cap in P1a and ∼ 5 pc from P4, meaning that illumination varies by a factor of ∼ 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 scoby5 .We link the spectral types of early-type NGC 6611 members published by Hillenbrand et al. (1993) to T eff and log g values using the tables of Martins et al. (2005), and then use the T eff and 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 ∼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 ≳1 pc away along the line of sight.We restrict our sample to stars within 2 pc in projection of the approximate cluster core (α, δ) J2000 = (18 h 18 m 35.s 9543, −13 • 45 ′ 20.364 ′′ ), which includes 30 stars with types B2.5 and earlier.We additionally restrict stars to have log 10 (L FUV /L ⊙ ) > 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 ∼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 ∼ 2500 Habing units towards the heads of P1a and P2, ∼ 1700 towards P3, ∼ 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 ∼ 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 ∼ 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 ∼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 ≳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ν > 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 ∼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 G0 We can roughly estimate the line-of-sight geometry of the Pillars using our two G 0 estimates, assuming that the catalog-based G 0, star is the emitted FUV radiation field and the FIR-based G 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 0, star and G 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 0, star = L/4π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 FIRbased estimate measures the apparent radiation field, so G 0, dust = L/4π(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.
(1) The heads of P1a and P2 are estimated to lie ∼0-1 pc and ∼2-3 pc, respectively, from the plane of the cluster; higher values are obtained when the background subtraction is applied to G 0, dust .The head of P3 is close to the size of the 160 µm beam, so beam dilution may slightly reduce G 0, dust and inflate separation; we estimate a separation of ∼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 ∼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 = Dcos(i) from the true distances D = 2.6, 3.2, and 3.6 pc and inclinations i = 47 • , 40 • , and 40 • 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 0, star ≥ G 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 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 ∼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 ∼50% brighter than the diffuse P4 emission in both bands and likely affects the G 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 0, dust is unsuitable for the LOS separation estimate.

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 pdrtpy6 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], 12 CO(J=1-0), 12 CO(J=3-2), and 12 CO(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 12 CO(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 12 CO(J=3-2), 2) there is zero emission outside the features and 3) the 12 CO(J=3-2) beam filling factor is the ratio of the beam areas Ω [C II] /Ω CO .We scale the 12 CO(J=3-2) intensities up by the ratio Ω CO /Ω [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 ≈ 2 × 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 ∼1-4 ×10 4 cm −3 .We estimate this uncertainty using a contour of a few times the minimum χ 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 ∼ 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 ∼ 1-2 ×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 selfconsistent across all modeled regions towards the Pillars.

Molecular and Atomic Hydrogen Column Densities
Molecular and atomic gas column densities are estimated using the 13 CO(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.

13 CO Column Density
We estimate the molecular gas column density assuming optically thin 13 CO(J=1-0) emission following Ti-   Mangum & Shirley (2015).We assume that the J=1-0 lines of 12 CO and 13 CO share the same excitation temperature along each line of sight, that 12 CO(J=1-0) is optically thick everywhere, and that the beam filling factor is unity.We adopt an isotopic ratio 12 CO/ 13 CO = 12 C/ 13 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 12 CO column density to molecular hydrogen column density N (H 2 ) using the abundance ratio 12 CO/H 2 = 8.5 × 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 Gas mass is calculated by summing N H over the boxes shown in Figure 12 and converting to mass using a mean molecular weight µ = 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 12 CO and 13 CO(J=1-0) observations and, for conversion of column density into mass, the heliocentric distance uncertainty, into the procedure described above.

C + Column Density
We do not detect the [ 13 C II] line towards the Pillars, so we estimate an upper limit on the C + column density using the [C II] channel noise T rms as a detection limit.We estimate an upper limit on the [ 13 C II] hyperfine component at +11.2 km s −1 with respect to the [ 12 C II] line (see Table 1 in Guevara et al. 2020) since it is sufficiently far from the [ 12 C II] line center given the observed linewidths, and we assume that the isotopic ratio 12 C + / 13 C + = 12 C/ 13 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 [ 12 C II] to [ 13 C II] given the [ 12 C II] optical depth.We find an upper limit of τ 12 ≲ 1.3 using the brightest [ 12 C 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 12 C + column density without a detected [ 13 C II] line described by Okada et al. (2015), and we assume a constant T ex rather than a constant optical depth.The adopted T ex = 107 K  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.
is calculated using the highest [C II] line brightness towards the Pillars and our upper limit of τ ≤ 1.3, and T ex is therefore a lower limit towards that location (considering the emitted radiation temperature T R constant, τ and T ex have an inverse relationship in the optically thin and τ ∼ 1 regimes).Assuming a fixed T 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 ex .T 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 ex ∼ 150 K, column densities would be approximately halved around the map, and it is unlikely that T ex is considerably lower than 100 K.We convert N (C + ) to atomic hydrogen column density N (H) using the abundance ratio C/H = 1.6×10 −4 (Sofia et al. 2004) and calculate column densities assuming a filling factor η = 1.We assume hydrogen to be entirely in the atomic phase where C + is present so that N H = N (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-σ 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 ∼30% of the molecular gas mass.

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 ≥250 µm maps have lower spatial resolution and do not resolve some Pillar features.We obtain ∼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 ext,160 /H.We take the R V = 3.1 value of C ext,160 /H = 1.9 × 10 −25 cm −2 /H from Draine (2003), following the method of Tiwari et al. (2021).This yields N 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 H and M H for 1000 realizations of the three uncertainties added to the observations and adopting the 1-σ uncertainties to be the difference between the 84 th and 16 th 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 H and M 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.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 AV is AV = NH/1.9× 10 21 cm −2 (Bohlin et al. 1978), so the column densities listed in the table can be divided by 1.9 to estimate AV .Column densities derived using the [C II] line have statistical uncertainty σstat = 0.6 (table units) except for location Shared-Base-E, which has σ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 nH and nH 2 are estimated using the C + and CO measurements of column density by dividing N (H) and N (H2), respectively, by the assumed LOS width of each feature.We assume NH = N (H) in the atomic gas and NH = 2N (H2) 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.
The dust-derived N 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 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.

Density
Line-of-sight averaged densities are derived by dividing the column density measurements by estimated lineof-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 ∼ 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 N 2 H + brightness temperatures (listed in Appendix A) according to estimates from the RADEX7 radiative transfer software (van der Tak et al. 2007).Densities ∼ 10 4 cm −3 are derived from C + -based column densities, which are consistent with those estimated using the PDR Toolbox.We adopt n H = 1.8 × 10 4 cm −3 as the atomic gas density by taking the average n 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 H2 = 1.3 × 10 5 cm −3 as the molecular gas density by taking the average n H2 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 H 2 ranges from 3-5 × 10 3 cm −3 at ∼100 K, so our adopted atomic gas density implies that the kinetic temperature T K should be similar to the excitation temperature T ex .Recalling that our T ex is a lower limit due to optical depth, we place a lower limit on atomic gas T K > 100 K.The derived molecular gas densities are greater than the critical density of 12 CO(J=1-0) by a factor of ∼ 50, so the radiation temperature T R , assumed equal to the observed T MB , of the optically thick 12 CO(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 12 CO must therefore be T K ∼ 50-70 K for most of the pillar bodies and reach T K ∼ 100 K towards brightly illuminated regions in P1a and P1b.The derived molecular gas density is close to the critical densities ∼ 10 5 cm −3 of HCO + and HCN and slightly below the critical densities ∼ 3 × 10 5 cm −3 of N 2 H + 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 (H 2 ), so their observed T MB = T R ∼ 15-20 K likely traces T K ∼ 20 K molecular gas towards P1a.

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 ionizationdriven shock has already been driven through the pillar heads.

Ionized Gas
Our observations do not probe directly the physical conditions of the ionized hydrogen (H II) phase, so we take values n e ∼ 1120-1800 cm −3 and T e ≈ 8000-10,000 K from optical ion line studies of the pillar surface by García-Rojas et al. (2006) andMcLeod 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 i /k B = 2n e T = 1.8-3.6×10 7K 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.

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 H2 = 1.3 × 10 5 cm −3 in Section 5.2.Molecular gas temperature should be coupled to dust temperature T d at densities n ≳ 10 5 cm −3 (White et al. 1999).We find T d ∼ 25 K towards the dense gas (Section 5.1.3),comparable to the value T ∼ 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 H2, therm /k B = n H2 T ≈ 3.2 × 10 6 K cm −3 .We estimate the pressure from turbulent support using the molecular line velocity dispersion σ obs ≈ 0.6 km s −1 (thermal velocity dispersion, σ therm ∼ 0.1 km s −1 , contributes less than 2% of σ obs ) to be P H2, turb /k B = ρσ 2 = (2n H2 µm H )σ 2 ≈ 1.5 × 10 7 K cm −3 (µ = 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 µ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 H2, B /k B ≈ 0.9-3.0 × 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 ⊥ ∼ 320 µG and P H2, B /k B ≈ 3.0 × 10 7 K cm −3 from our observed molecular gas velocity dispersion σ v ≈ 0.6 km s −1 along lines of sight with limited line-of-sight confusion.The magnetically-dominated total pressure within the molecular gas is P H2 /k B ≈ 4.8 × 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 ionizationdriven 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 H2, rot = ρω 2 r 2 where ωr ≈ 0.5 km s −1 based on the CO line velocities (see the PV diagram in Figure 9), finding P H2, rot /k B ≈ 10 7 K cm −3 (not included in our total pressure).This is comparable to turbulent pressure and less than magnetic pressure (3 × 10 7 K cm −3 , see above), and could contribute ∼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).

Atomic Gas
We determined the density of the warm, atomic gas to be n H ≈ 1.8 × 10 4 cm −3 .The lower limit on the T ex of the [C II] line indicates a kinetic temperature T 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 H /k B = n H T ∼ 2-5 × 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] σ therm ∼ 0.5 km s −1 given the assumed T K from the observed [C II] velocity dispersions σ obs ∼ 1.1-1.4km s −1 to find σ turb = (σ 2 obs − σ 2 therm ) 1/2 ∼ 1.0-1.3km s −1 .Adopting the turbulent velocity dispersion σ turb ≈ 1.3 yields a turbulent pressure of P H, turb /k B ∼ 5 × 10 6 K cm −3 in the atomic gas.This is an upper limit since the photoevaporative flow velocity ∼ 0.5 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 B ∼ 1 × 10 7 K cm −3 .The missing pressure is likely from the magnetic field in the atomic gas whose strength must be B ∼ 200µG.Pellegrini et al. (2007Pellegrini et al. ( , 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 ⊥ ∼ 85 µG and a corresponding P B /k B ∼ 2 × 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 τ AD ∼ 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 χ ion ∼ 10 −4 and τ AD ∼ 10 9 yr.

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 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 ⊙ , respectively, while the estimated H 2 masses are 72, 14, and 6 M ⊙ .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 where the degree of ionization χ ion is calculated for a typical cosmic ray ionization rate ζ CR = 3 × 10 −17 s −1 (Elmegreen 1979).For the adopted n H2 = 1.3 × 10 5 cm −3 , we find χ ion ∼ 6 × 10 −8 and τ AD ∼ 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 χ 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 τ AD ∼ 0.2 Myr, Note-Pressures in the molecular, atomic, and ionized gas near P1a.Pressures are expressed as P/kB 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.
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 massto-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.

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 ∼ 0.5 × 10 6 yr to photoevaporate the pillar, comparable to the ∼ 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 ∼ 10 5 yr on the pillar age required for unshocked pillar heads (White et al. 1999;Williams et al. 2001).Our estimated photoevaporation timescales ∼ 10 6 yr (Section 6.2) and the NGC 6611 cluster age ∼ 2 × 10 6 yr also favor the post-shock, ∼10 6 -yr-old Pillars scenario.See also the discussion of timescale tension by Williams et al. (2001).

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 ∼0.5 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 ∼0.1 Myr while the cloud is collapsing into a globule, and a wing during the stable cometary phase at an age of ∼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 12 CO 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 selfgravitating 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 12 CO 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.

Pillar Lifetimes
We estimate the photoevaporative lifetime of the pillars based on our observations.We adopt the equation for mass loss rate Ṁ = Av f ρ H from Gorti & Hollenbach (2002), where A is the total evaporative surface area, v f is the photoevaporative flow velocity from the surface, and ρ H is the mass density at the base of the flow.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 H = (P H,therm /ρ H ) 1/2 ∼ 1 km s −1 .The observed [C II] velocity dispersion gives an upper limit v f ≲ σ obs ≈ 1.4 km s −1 , though we attribute most of the observed velocity dispersion to turbulence (Section 5.3).We adopt v f ∼ 0.5 km s −1 .Mass density at the base of the flow ρ H = µm H n H , where we use µ = 1.33 and the atomic gas density n H = 1.8 × 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 N 2 H + 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 ∼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 ∼1-2 Myr by our estimate, consistent with the estimate by McLeod et al. (2015) of ∼3 Myr for these pillars as well as the estimate by Westmoquette et al. (2013) of ∼2 Myr for the pillars surrounding NGC 3603.This is similar to the ∼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 ∼Myr down to their present-day masses.As we discussed in Section 5.3, the photoevaporative lifetime is at least as large as τ 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 ∼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 ∼10.We do observe thinner (by ≲0.5 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 ∼60% from the Head to the Neck of P1a, and by ∼30% in P2.
We briefly consider the possibility that the necks do evaporate ∼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 ∼2.The meaningful difference between the two evolutionary scenarios is whether the Pillars are observed as globules or coherent pillars for the last ∼1 Myr of their lives; the work of Schneider et al. (2016) suggests the former.

CONCLUSION
We observed the Pillars of Creation in M16 in velocityresolved [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.
The largest pillar, P1, is the site of the brightest sion from both molecular gas and PDR tracers.It is separated into two substructures, P1a and P1b, which roughly correspond to its head and its base, respectively.These substructures are somewhat separated along the line of sight, though still linked by kinematically continuous [C II] emission.In P1a, the molecular tracers CO, HCN, HCO + , and CS reveal two dense filamentary structures which we term the "Threads" running parallel to the pillar body and terminating in a dense "Cap" at the head of the pillar.Our kinematic modeling indicates that these three distinct components are connected towards a region of high column density in the pillar head.Below the Threads past a region of low column density lies P1b, which is linked to the base of P2 by a bright, high column density PDR called the Shared Base.P2, a ∼1 pc long pillar which lies closer to the observer than P1, features a dense clump at its head and another halfway down its body.P2 and P3 are dwarfed in brightness by the PDRs towards P1.We include in our studies a handful of nearby PDR host structures, namely the Ridge and P4, which are kinematically related to the Pillars.
This wealth of data allows us to develop a geometric model for the structure of this pillar system, summarized by the schematic in Figure 10.Velocity gradients along the Pillars indicate that P1a and P3 point towards the observer while P2 points away.The relative strength of diffuse foreground ionized gas emission places P1a farthest from the observer, and P2 and P3 on the near side with respect to the bright O stars.P2 and P3 are darker optical sources and so their illuminated sides must face away from the observer.All of this is consistent with the conclusions of Pound (1998) andMcLeod et al. (2015).We establish from kinematic and column density information that P2 and P1b are connected by the Shared Base, a PDR which is significantly extended along the line of sight.
The [C II] line traces the same parsec-scale spatial and dynamic patterns as the molecular lines, exhibiting the same velocity gradients along P1 and P2 (P3 is too dim in [C II] and not well enough resolved to confirm or rule out a gradient).At smaller scales, [C II] and [O I] emission is generally smoother and more continuous than molecular emission, particularly towards the Threads in P1a, so we conclude that a PDR layer can manifest as a common envelope around dense features like the Threads.The [C II] and 12 CO(J=1-0) lines are blueshifted by ∼0.5 km s −1 w.r.t. the rest of the molecular lines towards the Shared Base and the heads of P2 and P3.We suggest that this is caused by the same process as the observationally established acceleration of gas away from the stars down the bodies of pillars (Pound 1998;McLeod et al. 2015) and that the process acts more strongly on the less dense gas in the outer layers of pillars.These motions are similar to those simulated by Lefloch & Lazareff (1994) and Mackey & Lim (2010) and may be a shear flow of the outer layers down the surface of the pillars.
We study the illumination based on both cataloged cluster members and FIR dust emission and constrain the FUV intensity to G 0 ∼ 10 3 Habing units towards the Pillars.The head of Pillar 1 is characterized by a relatively high G 0 ∼ 2000 Habing units and good agreement between the stellar and FIR estimates, indicating that it is directly illuminated and the projected distance to the stars is close to the true distance.
We determine the average molecular gas density in the Pillars to be ∼ 1.3 × 10 5 cm −3 and the average atomic gas density to be ∼ 1.8 × 10 4 cm −3 .The integrated column densities towards the Pillars yield total masses of 103, 78, 103, and 18 M ⊙ for P1a, P1b, P2, and P3 respectively.The atomic gas is typically ∼25% of the total mass, with the exception of P1b which contains about half its mass in the atomic phase.The masses of the Pillars are concentrated in the heads and occasional clumps along the bodies, such as the Horns atop P1b and the clump halfway down P2.These clumps are illuminated, as we detect PDR emission from their cluster-facing surfaces.
Based on pressures derived from our number densities, we assume that the pillars are in pressure equilib-rium with the ionized gas around them and that the PDRs associated with the Pillars are magnetically supported.In addition, the magnetic pressure of the dense molecular gas seems to be supporting the Pillar heads and clumps against gravity, but we infer that ambipolar diffusion is reducing the magnetic field support over a ∼1 Myr timescale.
We estimate that the heads of each pillar will be photoevaporated by radiation from NGC 6611 in ∼1-2 Myr.Whether they remain coherent pillars or evolve into free-floating globules as they evaporate may indicate whether or not flows of gas from the head into the body are a significant mass loss pathway compared to photoevaporation.Improved understanding of the evolutionary phases of pillars will reveal their relationships to cometary globules and other categories of observed features.
The Pillars' column shapes are to some extent a shadow of the concentrated mass at their heads, but it is clear from our observations and analysis that there are pre-existing density enhancements all along their bodies.The structure and orientation of the Pillars are affected by the pre-existing structure in the molecular cloud which is sculpted by the radiation field when massive stars turn on.3 and 8.The [C II] line spectra are shown with a dotted line before and a solid line after subtracting out the background as described in Appendix B; the second column in Table 3 lists whether the northern or southern background was subtracted.Vertical lines mark every 1 km s −1 between 20-28 km s −1 .Several lines/locations show additional components separated from the peak emission, such as CO at VLSR ∼ 29 km s −1 towards multiple positions or [C II] at VLSR ∼ 20-23 km s −1 towards the Ridge; all such components originate from background features identified in channel maps.Note-Peak main beam temperatures in Kelvins from lines at their native resolutions, which are listed alongside RMS temperatures in Table 1.The positions are shown in the leftmost panel of Figure 12.Since these are calculated as the maximum value from each spectrum, a measurement 1 or 2 times the RMS noise in Table 1 should be regarded as a potential nondetection; the spectra in Figures 13 and 14 provide additional context.The N2H + measurements are from the satellite line (J, F1, F) = (1-0, 0-1, 1-2) since the brightest lines at the center of the bandpass overlap (see Section 2).The Ridge is outside the half-power level of BIMA's primary beam, so the CO(J=1-0) measurements there are not given.The listed [C II] measurements have not had any background subtracted from them.The peaks of the background line profiles are thin (∼1-3 km s −1 ), which indicates that most of this emission originates from the atomic gas (Cuadrado et al. 2019).The ionized H II phase contributes ∼10% of the total [C II] emission on large (∼500 pc) scales (Tarantino On a more local scale around the Pillars, this diffuse emission is even brighter towards the southern end of P1, near P1b.Channel maps between V LSR ∼ 24-28 km s −1 indicate structure immediately east of P1b, as it is not clear if the structure is related to the Pillar system.This feature has a wider line profile than the northern background samples, as we see in the comparison in Figure 15.
In order to mitigate adverse effects of this background on our analysis of the Pillars, we conduct a background sampling and subtraction in order to help isolate [C II] emission from the Pillars themselves.We use the average of all pixels in four northern regions, drawn with the bold, solid line on the right panel of Figure 15, to correct "northern" targets such as P1a, P2, and P3.We use the average from within the southern region to correct spectra towards P1b and the Ridge; see the Background column of Table 7.These corrections are made for all [C II] spectra we show going forward and particularly throughout the kinematic analysis in Section C where the shape of the line is important to our analysis.We do not account for optical depth while subtracting background spectra; we find in Section 5.1.2that [C II] has an optical depth ≲ 1, so adverse effects of background subtraction on line shape should be minimal.We do not make a channel-by-channel spectral subtraction while calculating C + column densities in Section 5.1.2.Some distinct components are identified in the background spectra from Regions #2 and #3 around 30-40 km/s in Figure 15.Referencing the channel maps in Figure 4, these can be associated with a morphologically distinct north-to-south strip of redshifted gas which does not appear directly associated with the Pillars.The influence of this redshifted component is not seen in the averaged background spectrum, and the feature itself has very little overlap with P1, so it should not adversely affect our background subtraction.
The background subtraction proves to be useful in determining whether characteristics of the Jackknife tests using combinations of these and other background samples confirm that the choice of background samples doesn't upset results, i.e. change the profile and velocity of lines too much, but reveal limitations in our ability to make specific and precise claims about background-subtracted [C II] line intensities since the background varies throughout the system (see the velocity shift between Regions #1 and #2 in Figure 15).

C. KINEMATICS AND GEOMETRY OF P1A
P1a is composed of three morphologically and kinematically distinct components, making it the most complex structure in the pillar system.A 3D geometrical model, informed by an in-depth kinematic analysis of P1a, is necessary to contextualize other physical characteristics derived from our observations.We fit Gaussian line profiles to observed [C II] and molecular line spectra in order to quantify, or at least approximate, characteristics of the observed profiles.In advance of this kinematic modeling, we subtract channel-by-channel from [C II] spectra the background identified near the Pillars.We describe the identification and correction in Appendix B. Since we only model spectra towards P1a, we only use the northern background sample described in Appendix B.

C.1. Kinematic Modeling Towards P1a
We fit 1, 2, 3, and 4 component models to each HCO + and [C II] spectrum in a ∼ 1' box surrounding P1a (∼500 [C II] spectra, ∼10,000 HCO + spectra).For each line, all spectra were modeled using the same model template and initial conditions.The strength of this unsupervised pixel grid fit is that real, physically meaningful patterns may emerge from this large number of modeled spectra.We visualize the fitted models in (RA, Dec, Velocity) space by plotting each pixel's fitted component line centers.There are 1-4 fitted components per pixel depending on the model template and the number of components with sufficiently large amplitudes with respect to the noise.We plot these components using a 2D histogram by projecting along Dec to view the RA-Velocity plane; due to the components' spatial orientations, we find the RA-Velocity projection easier to interpret than the Dec-Velocity projection and so we limit our discussion to the former.These histograms are similar in orientation and meaning to PV diagrams.
These projections, shown in Figure 16, all illustrate a spread of components towards the east which move kinematically closer together towards the Merge Point in the west, where they appear to merge together in the [C II] series and the HCO + 1-component projection.In the HCO + model series, three line center groupings are observed towards the eastern side of the image: one at V LSR ≈ 25.5 km s −1 , one at 24.5-25 km s −1 , and one at 23.5 km s −1 in the east with a strong gradient towards higher velocity to the west.These velocities and their gradients are consistent with those of the two Threads and the Cap which we observed in the channel maps of most lines.These groupings are observed in the 1, 2, and 3 component figures, and these groupings do not fundamentally change or disappear even when a 4th component is made available.In the [C II] model series, we only observe two groupings at V LSR ≈ 23.5 km s −1 and 25.5 km s −1 , and as with HCO + this pattern persists even when more than two components are allowed.
The HCO + results show a shift in behavior towards the western part of the head, near the Merge Point in Figure 1, between 1 and 2 allowed components.In the 1-component models, we observe the kinematic merging of the high-and low-velocity component groupings into a single group around V LSR ≈ 24.7 km s −1 .When we make a second component available to the HCO + line models, we see the three eastern velocity groups converge towards the same velocity but then dramatically shift into a two-component grouping towards the Merge Point where the 1-component models converged.This pattern persists even as 3 and 4 components are available to the HCO + spectra.This is peculiar for two reasons: first, it is at odds with the picture of several components merging together which we see in the HCO + single component model and all the [C II] models; and second, it relies on an abrupt shift in component velocity coupled with the abrupt disappearance of an entire component.We suspect that, rather than a physical phenomenon, this is a sort of computational "phase shift" within the solution space of the model fit in which the result abruptly transitions from one solution to another between nearby pixels.A more detailed investigation of the affected spectra is required in order to tell if the components really do merge and the HCO + "phase shift" is really just a computational artifact, or if the situation too complex for us to make such a claim.For a few locations towards the head, we perturbed the models' initial conditions and examined the fitted models to check whether our results are robust and whether we come to the same conclusions as described above.Through this more supervised method, we find that all results from the unsupervised method are sound except those towards the northwestern corner of the head around the Merge Point, where the solution space phase shift occurs.We find that the 2-component solutions towards that location are not unique and that a single component works comparably well.We conclude that HCO + line spectra towards the eastern and southern parts of the head can be decomposed into three distinct components in all the molecular lines, and the components are consistent from line to line.The [C II] spectrum towards these locations can be decomposed into two components which generally correspond to the highest and lowest velocity molecular gas components (the Eastern Thread and the Cap).HCO + line spectra towards the Merge Point can be decomposed into 1-3 components with similar reduced χ 2 and so a single component is the least complex solution.The three components merge into one towards the Merge Point, as we see in the single-allowed-component panel in Figure 16, and become too indistinct to fit separately.

C.2. Geometry of P1a
We interpret these results to originate from a pillar head composed of three distinct molecular gas components, the two Threads and the Cap, embedded in a warm PDR gas envelope.We observe the two Threads as spatially separate entities south of the pillar head and we can distinguish the components in the eastern and southern regions of the head in velocity, even though they are spatially blended.Due to its strong velocity gradient, the Cap is the most kinematically separated from the rest of the head to the east.At that location, we observe the widest line profiles which we are able to decompose into three components.These spectra can also be fit with two components, where one represents the Cap and the other wider component represents a merged-Thread component, but they cannot be well fit with a single component.
These three morphologically distinct components spatially overlap towards the head and kinematically overlap towards the northwest part of the head, indicating that they are physically merged towards the northwest.The gradients of the components are similar in absolute value but not in sign and are organized radially around the Merge Point.Line profiles in every line towards the Merge Point are consistent with a single emitting component, while line profiles towards other locations of the head cannot be modeled by a single component.We detect N 2 H + and C 18 O emission and high 13 CO column density towards the Merge Point, and all lines peak in brightness near that region.The N 2 H + detection implies sufficiently cold, shielded gas in addition to high column density, which is more likely towards a physical intersection of clouds rather than a projected stack of clouds.
Component interactions and overlays like this one may be responsible for the broad [C II] line profiles observed towards other regions.Ossenkopf et al. (2013) describe "macro-turbulent" motions, which cause a greater portion of the [C II] line to be optically thick and create a broader line profile than expected for opacitybroadening.Our line profile analysis and high spatial resolution observations reveal the individual com-ponents towards P1a, but they are likely obscured by geometry or insufficient spatial resolution elsewhere.
The [C II] emission towards P1a shares broad kinematic characteristics like the gradients along the Threads with the molecular line emission as we see in Figure 16, but some key differences noted in Section 3.2 lead us to conclude that the PDR gas forms an extended envelope around the molecular gas.We do not detect the threaded morphology just south of the pillar head in [C II] or [O I] (see Figure 3).This is not solely due to the coarser spatial resolution of the fine structure line observations because smoothing the CO to the same resolution preserves the structure.The [C II] emission along the pillar body towards the Threads is spatially broad, ∼0.3 pc, and is centered on the Eastern Thread; from a purely spatial standpoint, the [C II] emission appears to trace only the Eastern Thread.We do not detect the same transverse velocity gradient across the Eastern Thread in [C II] as we do in the molecular lines.The PV diagrams in Figure 9 show that [C II] shifts to lower velocities instead, such that the [C II] line peaks at similar velocities as the molecular lines towards both the Eastern and Western Threads.We interpret this as a kinematic detection of [C II] emission associated with the Western Thread.Modeled [C II] central velocities show that, while the Eastern Thread has a bulk molecular gas velocity ∼1 km s −1 higher than the Western Thread, the [C II] emission only shifts by ∼0.3 km s −1 between the same locations.The [C II] PV diagrams are spatially and kinematically broad and uniform, and the molecular line PV diagrams show more structure with more extreme velocities and narrower line widths.

Figure 1 .
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 ′ ≈ 0.5 pc at 1740 pc line of sight.
3. RESULTSWe present a rich collection of observations of the Pillars of Creation alongside publicly available archival data.Integrated intensity and continuum maps in Figure2show 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.
The [C II] and [O I] line profiles (FWHM ∼ 2-3 km s −1 ) match the expected line width for atomic gas PDR emission.The peak main beam temperatures T MB of [C II] reach ∼40 K towards the Pillars.The 12 CO line widths ∼ 2 km s −1 and peak T MB (reaching ∼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 ∼ 1-2 km s −1 ) and reach lower T MB (∼20 K in HCO + and HCN; ∼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).

Figure 2 .
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 VLSR = 18-27 km s −1 , as they tend to have longer blue wings, and molecular lines not otherwise specified are integrated between VLSR = 19-27 km s −1 .N2H + and HCN are integrated between VLSR = 12.6-32 km s −1 with respect to the rest frequencies given in Table1in order to include several satellite lines for each species (see Section 2).Contours are [0.5, 3.0] K km s −1 for C 18 O(J=1-0); [4, 24, 44] K km s −1 for 13 CO(J=3-2); and [9, 18, 27, 45] ×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 .
Figure 3. Integrated [C II] (black) and CS (white) line intensities between VLSR = 19-27 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 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 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 largescale [C II] emission indicated by the lowest two contours is the VLSR ∼ 25 km s −1 background discussed in detail in Appendix B.
Figure 3. Integrated [C II] (black) and CS (white) line intensities between VLSR = 19-27 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 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 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 largescale [C II] emission indicated by the lowest two contours is the VLSR ∼ 25 km s −1 background discussed in detail in Appendix B.

Figure 4 .
Figure 4. 158 µm [C II] line channel maps binned to 1 km s −1 and centered on the Pillars.The beam is shown in the lower right corner of each map.

Figure 6 .
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 km s −1 between 20-27 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 km s −1 .The [O I] contours are placed at [0, 6, 12, 18, 24, 30] K 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).

Figure 7 .
Figure 7. HCO + (color) and CS (contour) line channel maps binned to 1 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 .
Figure 8. Position-velocity (PV) diagram along each pillar using the [C II] and 12 CO(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 TMB in Kelvins.Emission at VLSR ≳ 27 km s −1 is from background features not directly related to the Pillars.Velocities are VLSR.

Figure 9 .
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.

Figure 11 .
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 χ 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.

Figure 12 .
Figure 12.The top-left, top-right, and bottom-left panels show hydrogen column densities NH 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 Table4, and red boxes show the same for Table6.Orange squares mark the locations from which column densities are sampled in Table3.In Table3, 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.

Figure 13 .
Figure 13.Spectra observed towards the pillars in several lines.All data are convolved to the [C II] resolution, except CO(J=3-2) which is shown at its native resolution.Locations are those shown in the leftmost panel of Figure 12 and listed by name in Tables3 and 8.The [C II] line spectra are shown with a dotted line before and a solid line after subtracting out the background as described in Appendix B; the second column in Table3lists whether the northern or southern background was subtracted.Vertical lines mark every 1 km s −1 between 20-28 km s −1 .Several lines/locations show additional components separated from the peak emission, such as CO at VLSR ∼ 29 km s −1 towards multiple positions or [C II] at VLSR ∼ 20-23 km s −1 towards the Ridge; all such components originate from background features identified in channel maps.

Figure 14 .
Figure 14.Same as Figure 13 showing different species.All data except 13 CO(J=3-2) are convolved to the [C II] resolution.The line intensity axis spans a smaller range here.Multiple satellite lines appear in the HCN spectra.The N2H + line is the satellite line (J, F1, F) = (1-0, 0-1, 1-2).The CO emission around VLSR ∼ 29 km s −1 is from a background feature.kept at least one beam (15 ′′ ) away from the optical/NIR edges of the Pillars, which ensures that the regions don't overlap with the Pillars and include Pillar emission while remaining close to the Pillars so that the background spectra are relevant to the Pillars.The [C II] spectra from within these regions include significant emission

Figure 15 .
Figure 15.(Left) The [C II] integrated intensity between VLSR = 18-27 km s −1 , the same interval shown earlier in Figure 2, is shown in color.The four northern background regions (numbered) and the single southern background region are outlined.(Right) The mean spectrum from within each individual background region on the left panel, labeled accordingly.The "Average" background spectrum drawn in the bold, solid line is the average across all pixels in the numbered northern regions and does not include the southern region.Vertical lines mark every 1 km s −1 between 20-30 km s −1 .et al. 2021).[C II] line profiles from the diffuse, ambient H II region should have ∼10 km s −1 widths reflective of the turbulent velocity dispersion observed by Higgs et al. (1979) towards the M16 H II region in lines of H and He.The background [C II] line profiles have wings that could arise from a wide, low-intensity component originating from the ambient H II region.This will be explored in greater detail in a future study.The characteristic PDR velocity V LSR ∼ 25-26 km s −1 of the region surrounding the Pillars may originate from precluster dynamics and poses a unique challenge for the analysis of P1, which lies almost exactly at this velocity.On a more local scale around the Pillars, this diffuse emission is even brighter towards the southern end of P1, near P1b.Channel maps between V LSR ∼ 24-28 km s −1 indicate structure immediately east of P1b, as it is not clear if the structure is related to the Pillar system.This feature has a wider line profile than the northern background samples, as we see in the comparison in Figure15.In order to mitigate adverse effects of this background on our analysis of the Pillars, we conduct a background sampling and subtraction in order to help isolate [C II] emission from the Pillars themselves.We use the average of all pixels in four northern regions, drawn with the bold, solid line on the right panel of Figure15, to correct "northern" targets such as P1a, P2, and P3.We use the average from within the southern region to correct spectra towards P1b and the Ridge; see the Background column of Table7.These corrections are made for all [C II] spectra we show going forward and particularly throughout the kinematic analysis in Section C where the shape of the line is important to our analysis.We do not account for optical depth while subtracting background spectra; we find in Section 5.1.2that [C II] has an optical depth ≲ 1, so adverse effects of background [C II] line profile should be associated with the Pillars.The molecular and (unsubtracted) [C II] spectra through the pillar bodies, and particularly through the head of P1 in Figure 13, differ noticeably in that the [C II] line profile has more emission at higher velocities than the other tracers.No shift is detected in the [C II] line center w.r.t. the molecular lines, so the asymmetry in the line profile must be limited to the low-intensity line wings.If we were to associate all [C II] emission towards the head of P1 with the pillar itself, then we would conclude that the PDR layer has a significantly different kinematic signature than the deeper molecular gas layers, thus resulting in an extended red tail in the [C II] spectrum.However, when we apply the [C II] background subtraction, we find that the redshifted [C II] tail disappears and the [C II] spectrum looks much more like the molecular lines, suggesting that the dynamics of the PDR and molecular layers are not dissimilar.The background identification and subtraction is therefore necessary in order to identify Pillar-related spectral characteristics and avoid ascribing background emission to the Pillars.

Figure 16 .
Figure 16.RA-Velocity projections of the 1, 2, and 3 component line center solutions.4 component model results are not shown here but are not significantly different than the 3 component results.

Table 2 .
Fitted G0 and n to observed line ratios using pdrtpy.

Table 3 .
Column and number densities.

Table 4 .
Atomic and molecular gas mass.Note-Gas mass M derived for the pillars from [C II],13 CO(J=1-0), and dust emission.Masses are given in M⊙ 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.

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

Table 6 .
Pillar head and neck masses, mass loss rates, and photoevaporative lifetimes.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 nH = 1.8 × 10 4 cm −3 .Adopted flow velocity is 0.5 km s −1 .

Table 8 .
Peak TMB for each line.