<figure><img src="image_1.png" /></figure>

# Gas Kinematics and Dynamics of Carina Pillars: A Case Study of G287.76-0.87

Ngo-Duy Tung<fnref n="1" />,<fnref n="2" />,<fnref n="3" /> ID, Le Ngoc Tram<fnref n="4" /> ID, Archana Soam<fnref n="5" />,<fnref n="6" /> ID, William T. Reach<fnref n="6" />,<fnref n="7" /> ID, Edwin Das<fnref n="8" />, Ed Chambers<fnref n="6" />,<fnref n="7" /> ID, Blesson Mathew<fnref n="8" /> ID, and Heiko Richter<fnref n="9" /> ID

## Abstract

We study the kinematics of a pillar, namely G287.76-0.87, using three rotational lines of ¹²CO(5-4), ¹²CO(8-7), ¹²CO(11-10), and a fine structure line of [O I] 63 µm in southern Carina observed by SOFIA/GREAT. This pillar is irradiated by the associated massive star cluster Trumpler 16, which includes η Carina. Our analysis shows that the relative velocity of the pillar with respect to this ionization source is small, ~1 km s⁻<fnref n="1" />, and the gas motion in the tail is more turbulent than in the head. We also performed analytical calculations to estimate the gas column density in local thermal equilibrium (LTE) conditions, which yields N<sub>CO</sub> as (~0.2–5) × 10¹<fnref n="7" /> cm⁻<fnref n="2" />. We further constrain the gas’s physical properties in non-LTE conditions using RADEX. The non-LTE estimations result in n<sub>H₂</sub> ≃ 10⁵ cm⁻<fnref n="3" /> and N<sub>CO</sub> ≃ 10¹<fnref n="6" /> cm⁻<fnref n="2" />. We found that the thermal pressure within the G287.76-0.87 pillar is sufficiently high to make it stable for the surrounding hot gas and radiation feedback if the winds are not active. While they are active, stellar winds from the clustered stars sculpt the surrounding molecular cloud into pillars within the giant bubble around η Carina.

*Unified Astronomy Thesaurus concepts: Interstellar clouds (834); Interstellar medium (847); Photodissociation regions (1223)*

## 1. Introduction

The Carina Nebula, one of the largest and brightest diffuse nebulae in the southern sky, is home to many massive star-forming processes. In this region of star formation, the most massive stars (O/B stars) generate intense ultraviolet (UV) radiation of energy E > 13.6 eV as feedback to the nascent molecular cloud (see Figure 1). These energetic UV photons ionize atomic hydrogen (H) and create an “H II region” (see, e.g., Strömgren 1939; Kahn 1954) around the stars. At the boundaries of the H II region, the UV photons that are less energetic than the ionization potential of H (that is, E < 13.6 eV) escape and penetrate the surrounding molecular gas. Along the way, they dissociate molecules, heat dust grains, and ionize atoms such as carbon, silicon, and sulfur, which have ionization potentials lower than 13.6 eV. This gives rise to a region where H is neutral, but low-ionization-potential atoms are ionized, and molecules such as CO are photodissociated, which is called a photodissociation region or PDR. (Hollenbach & Tielens 1997).

At the interface between a molecular cloud and an H II region, there exist some of the most spectacular features, namely globules and pillars. Pillars, which arise when UV radiation from hot stars propagates into a dense, clumpy region in the surrounding molecular cloud (Hester et al. 1996), have an elongated column-like shape and a physical connection to the gas reservoir of the cloud. Globules, which result from the separation

of pillars from the gas (White & Gee 1986; White 1993; Thompson et al. 2004; Schneider et al. 2012), are isolated (often with a young star located near the head) and have a head–tail structure pointing toward the illuminating source. In pillars and globules, the gas is strongly influenced by radiation; thus, the absorption of UV photons of E < 13.6 eV by the gas leads to intense emission of [C II] at 158 µm and [O I] at 63 and 145 µm and H₂ rovibrational transitions. Deeper in PDRs, CO rotational and [C I] 370 and 609 µm lines can be observed. Therefore, observations of these spectral lines in high spectral resolution (Mookerjea et al. 2019; Anderson et al. 2019; Kirsanova et al. 2020) can help us to trace the gas kinematics and understand the physical conditions of the gas in the emitting region.

The energy injection from an intense ionization source plays both positive and negative roles in the surrounding clouds. This feedback could trigger a new generation of stars to form (Elmegreen & Lada 1977) or halt this process or even disrupt the cloud structure (Pabst et al. 2019). In the case of the southern Carina, Rathborne et al. (2004) showed that a cluster of young stellar objects are likely to form in the pillar G287.84-0.82 (known as Treasure Chest; see Figure 1), but there is no evidence of star formation in the neighboring pillar (G287.76-0.87, highlighted by a box in Figure 1).

In this work, we will focus on the latter pillar, which is shown in Figure 2. We aim at characterizing the kinematic structure of this cloud using the very high spectral resolution of mid-infrared rotational transitions of CO and the fine-structure line of oxygen observed by the German Receiver for Astronomy at Terahertz Frequencies (GREAT) instrument (Risacher et al. 2018) on board the Stratospheric Observatory for Infrared Astronomy

<figure><img src="image_2.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ University of Science and Technology of Hanoi, Vietnam Academy of Science and Technology, 18 Hoang Quoc Viet, Hanoi, Vietnam</footnote>

<footnote n="2">² Université de Strasbourg, CNRS, Observatoire astronomique de Strasbourg, UMR 7550, 11 rue de l'Université, 67000, Strasbourg, France</footnote>

<footnote n="3">³ Université Paris Cité, Université Paris-Saclay, CEA, CNRS, AIM, F-91191, Gif-sur-Yvette, France</footnote>

<footnote n="4">⁴ Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, 53-121, Bonn, Germany; nle@mpifr-bonn.mpg.de</footnote>

<footnote n="5">⁵ Indian Institute of Astrophysics, II Block, Koramangala, Bengaluru 560034, India; archana.soam@iiap.res.in</footnote>

<footnote n="6">⁶ SOFIA Science Center, Universities Space Research Association, NASA Ames Research Center, Moffett Field, CA 94035, USA</footnote>

<footnote n="7">⁷ Space Science Institute, 4765 Walnut Street, STE B, Boulder, CO 80301, USA</footnote>

<footnote n="8">⁸ Department of Physics and Electronics, CHRIST (Deemed to be University), Bangalore, Karnataka 560034, India</footnote>

<footnote n="9">⁹ German Aerospace Center (DLR), Institute of Optical Sensor Systems, Rutherfordstraße 2, 12489 Berlin, Germany</footnote>

<footnote>Received 2023 January 16; revised 2024 March 29; accepted 2024 March 31; published 2024 July 4</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Southern Carina Nebula observed by Spitzer at 8 µm. Our target (pillar G287.76-0.87) is highlighted by a green box and located next to the Treasure Chest pillar (G287.84-0.82, Smith et al. 2005). The ionization source is the Tr 16 cluster, including the supermassive star η Carina, shown by a star symbol. The area around η Carina is white because it saturates the 8 µm camera.</figcaption></figure>

(SOFIA; Young et al. 2012). Furthermore, being the cooling lines of PDRs, these lines are ideal tracers to constrain the pillar's physical properties.

This paper is structured as follows. Our study's observations and data acquisition are presented in Section 2. The results obtained and the analyses applied to the results are presented in Section 3. We discuss our results in Section 4. Finally, a summary and conclusion of our findings are given in Section 5.

## 2. Observations and Data Acquisition

The observations used in this work are made using the GREAT on board SOFIA, in its 4GREAT/High Frequency Array (HFA) configuration. GREAT is able to observe five frequencies simultaneously. The 4GREAT channel, which consists of four coaligned pixels at four frequencies (Durán et al. 2021), was tuned to ¹²CO(5-4), ¹²CO(8-7), ¹²CO(11-10), and ¹²CO(22-21). The HFA, which consists of seven independent pixels in a hexagonal geometry with a central beam, was tuned to the [O I] line at 63 µm at 4.7 THz.

Observations of G287.76-0.87 were made during the 2019 SOFIA deployment to New Zealand as part of the Director's Discretionary Time project (75_0038). The data were obtained in about 130 minutes, spread across two flights on 2019 June 26 and 27. The maps were obtained in single beam-switch mode and designed to fully sample the 6.3 [O I] beam, with

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Zoom-in on the pillar G287.76-0.87: the dashed lines indicate the cuts for PV diagrams presented in Section 3.3, and the white circles mark their corresponding zero positions.</figcaption></figure>

spectral readouts every 3". Each readout had 0.4 s of on-source integration time. For the [O I] map, each of the 7 pixels creates its own independent map, resulting in about 2.8 s of on-source time per map point. For the 4GREAT maps, the 3" spacing resulted in oversampled maps. During the observations, the 4G4 channel that was tuned to CO(22-21) was unstable, and therefore, those data are not included in our analysis.

The data were reduced and calibrated by the GREAT team using their kalibrate software package (Guan et al. 2012) for atmospheric calibration and GILDAS/CLASS<fnref n="10" /> for final beam calibration, baseline removal, and map convolution and gridding. The final maps have pixel sizes of 18", 8.6", 6.6", and 2" for the CO(5-4), CO(8-7), CO(11-10), and [O I] maps, respectively. The spectra were all smoothed from their native resolution to channel widths of ~0.3 km s⁻¹. All data are publicly available at the Infrared Science Archive.<fnref n="11" />

## 3. Results

### 3.1. Channel Maps

We investigated the motion of gas in the pillars using velocity channel maps. Figure 3 shows the channel maps of [O I] 63 µm, whose spatial resolution is the highest in our data set. The channel of the CO lines will be shown later in Figures 11 and 12 in Appendix A. As the velocity increases, the signal becomes stronger, peaks at −14.7 km s⁻¹, and then starts to fade. The velocity for all four lines is negative, which means that the gas motion is blueshifted with respect to the local standard of rest (LSR) by ~−14 km s⁻¹. This velocity is redshifted with respect to the nearby cloud in η Carina that has the LSR velocity ~−20 km s⁻¹. We discuss the motion of the pillar in Section 4.1. The southwest region of the gas seems very peculiar in [O I] 63 µm, which is not seen in CO transitions (Figures 11 and 12 in Appendix A) owing to the larger beam sizes of the CO observations.

<footnote n="10">¹⁰ http://www.iram.fr/IRAMFR/GILDAS</footnote>

<footnote n="11">¹¹ https://www.sofia.usra.edu/node/2684</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Channel maps of [O I] 63 µm, the highest spatial resolution in our data set. The motion is blueshifted, indicating that the pillar is in front of the radiation source ($\eta$ Carina).</figcaption></figure>

## 3.2. Moment Maps

The upper left panel of Figure 4 shows the integrated brightness map of the [O I] 63 µm line. The map shows that the emission peaks at the northwest of the pillar and decreases radially along the northwest-to-southeast direction. This can be reconciled with the fact that the northwest region of the pillar is exposed to the radiation source ($\eta$ Carina). Radiation and interaction with the expanding H II shell can excite the molecular and atomic gas, and excitation becomes radially weaker behind this interaction front.

Radiation and the expanding shell H II also heat dust grains. Therefore, we compare the spatial correlation between gas emission and dust continuum in the pillar G287.76-0.87. The results are shown in Figure 4 (upper right and lower panels). The contours are the integrated intensity of [O I] (yellow) and CO (other colors) emissions. The background represents the 8 µm Spitzer map, which can trace the emission from polycyclic aromatic hydrocarbon (PAH) molecules. The [O I] 63 µm emission correlates very well with the 8 µm PAH continuum emission and is located ahead of the CO lines, representing the PDR structure.

The left panel of Figure 5 shows the first moment map of the [O I] 63 µm lines, in units of kilometers per second. The maps for the CO lines are shown in Figure 14 in Appendix B. To reduce noise, a mask is applied to these maps to select the data within the zeroth moment contours as shown in Figure 4. Along the line of sight, the gas velocity at the head is more negative than that at the tail of the pillar.

In the right panel of Figure 5, we show the second-moment map of [O I] 63 µm, in units of kilometers squared per second. Here we choose not to show the moment maps for the other CO tracers because of their low spatial resolutions, which result in little meaningful information. Interestingly, although

having a lower intensity-weighted velocity as suggested by the first-moment maps, the tail's motion is more turbulent than that of the head, with a maximum dispersion of around 1–3 km s⁻¹. A very similar effect has recently been noted in stellar wind shells, specifically the Orion Veil (Kavak et al. 2022) and RCW 120 (Kirsanova et al. 2023).

## 3.3. Position–Velocity Diagrams

With all information from the channel maps and moments maps, we can reconstruct a 3D motion of the pillar if we know how the velocity is related to the position. In this section, we aim at that goal by means of position–velocity (PV) diagrams with a few slices through the pillar. Using the Python module *pvextractor* (Ginsburg et al. 2016), Figure 6 shows the PV diagrams in CO(5-4), CO(8-7), CO(11-10), and [O I] 63 µm along the three cuts shown in Figure 2.

The first cut is set to point toward the $\eta$ Carina. The zero position for this cut is at (10$^h$45$^m$21$^s$, −59°58′44). Apparently, CO(5-4), CO(8-7), and [O I] 63 µm show a clear structure, whereas CO(11-10) does not because this transition does not trace well in this direction (see Figure 4, bottom right panel). Signals rise sharply in the northwest and narrow down in the southeast.

The second cut is perpendicular to the first. The zero position for this cut is at (10$^h$45$^m$22$^s$, −59°58′44). CO(5-4) seems to be compact with little velocity spread, and the CO(8-7) shows that the gas motion is not highly dispersed, yet the spatial spread of the gas at the same velocity is quite visible. On the other hand, CO(11-10) lines are broader at this position, with FWHM of 4 km s⁻¹. The gas motion spans a range of ~5 km s⁻¹, which is larger than those of other tracers. The plot for [O I] 63 µm in this cut does not provide much information. It seems to have more noise compared to the main signal.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Integrated intensity maps of [O I] is shown in the upper left panel. Other panels show the correlation between the [O I], CO(5-4), CO(8-7), CO(11-10) integrated intensities (contours) and the Spitzer 8 µm map (background color).</figcaption></figure>

The third cut, whose zero position is at (10$^{th}$45$^{m}$23$^{s}$, -59°59'75), is made based on the observation that there is a strange feature (a hole) at the southeast of the clump. In this cut, the velocity component is not seen by CO(11-10) and [O I] 63 µm, suggesting that the gas is more diffuse than the rest of the clump, as we saw in the Spitzer image (Figure 2).

### 3.4. Spectral Distribution

Figure 7 shows the average spectra for all the lines. The velocity resolution of 4GREAT is finer than 0.1 km s⁻¹ (Durán et al. 2021), which well resolves the spectral lines. We then fit these spectra with a single Gaussian function. The corresponding mean µ and standard deviation σ are marked in each panel. For each line, our confidence interval is selected between µ - 3σ and µ + 3σ, which accounts for 99.7% of the signal. For the CO(8-7) line, since a single Gaussian function does not seem to cover all the signal, we try to fit it with a double Gaussian function. The parameters for the second Gaussian are µ = -8.1 km s⁻¹ and σ = -0.97 km s⁻¹. The nature of this second component needs to be determined in future studies.

From these fitted parameters, one can see that the [O I] line is broader than the CO lines.

### 3.5. Physical Properties

Our goal in this section is to constrain some key physical properties of the pillar G287.76-0.87, i.e., the kinetic temperature T<sub>kin</sub>, the number density n<sub>H₂</sub>, and the column density N<sub>CO</sub> (which can convert to the gas column density), using the observed spectra. CO transitions may be either optically thick or optically thin, but we do not have data for other isotopes of CO (e.g., ¹³CO, which is optically thin) and/or C⁺ (e.g., ¹⁵C⁺), to measure optical depth. Therefore, we first propose a different analytical approach to estimate the gas column density in local thermal equilibrium (LTE) conditions and then perform the non-LTE modeling using the radiative transfer program RADEX.

#### 3.5.1. LTE Calculations

Our data set contains only three transitions of the ¹²CO and one atomic [O I] transition. The ¹²CO transitions can be either
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Left: first-moment map of [O I] 63 µm integrated in velocity according to the channel as in Figure 3. Along the LOS, the head velocity is more negative than the tail velocity. Right: corresponding second-moment map of [O I]. The motion projected on the LOS in the tail is shown to be more turbulent than in the head.</figcaption></figure>

optically thick or optically thin, whereas the [O I] 63 µm seems to be thick in the PDR (see Goldsmith 2019). In this section, we, therefore, estimate the physical properties from the line ratios of $^{12}$CO in LTE conditions by assuming that they are either optically thin or optically thick, as introduced in Hsieh et al. (2012), with the condition that all the lines are tracing the medium with the same excitation temperature, and their filling factors are all equal.

Let us consider a rotational transition at frequency $\nu$. Its brightness temperature in LTE (i.e., the source function is approximated as a blackbody function) is

$$T_B = f[J_{\nu}(T_{\text{ex}}) - J_{\nu}(T_{\text{bg}})][1 - \exp(-\tau_{\nu})], \quad (1)$$

where $f$ is the filling factor, $T_{\text{ex}}$ is the excitation temperature, $T_{\text{bg}} = 2.73$ K is the background temperature, $\tau$ is the optical depth, and $J_{\nu}$ is the blackbody temperature defined at temperature $T$ as

$$J_{\nu}(T) = \frac{h\nu/k_{\text{B}}}{\exp(h\nu/k_{\text{B}}T) - 1}. \quad (2)$$

In the first case where both lines are optically thick ($\tau_{\nu} \gg 1$), the ratio of the integrated intensity of line 1 ($\nu_1$, $\tau_{\nu_1}$) over that of line 2 ($\nu_2$, $\tau_{\nu_2}$) is

$$R_{\nu_1/\nu_2} = \frac{[J_{\nu_1}(T_{\text{ex}}) - J_{\nu_1}(T_{\text{bg}})]}{[J_{\nu_2}(T_{\text{ex}}) - J_{\nu_2}(T_{\text{bg}})]}. \quad (3)$$

In the second case, both lines are optically thin ($\tau_{\nu} \ll 1$), which yields

$$R_{\nu_1/\nu_2} = \frac{[J_{\nu_1}(T_{\text{ex}}) - J_{\nu_1}(T_{\text{bg}})]\tau_1}{[J_{\nu_2}(T_{\text{ex}}) - J_{\nu_2}(T_{\text{bg}})]\tau_2}, \quad (4)$$

with $\tau$ being the opacity, and the opacity ratio

$$\frac{\tau_1}{\tau_2} = \frac{\nu_2^2 g_{u1} A_1 N_1^{\text{tot}} Q_2 e^{-E_{l1}/T_{\text{ex}}}}{\nu_1^2 g_{u2} A_2 N_2^{\text{tot}} Q_1 e^{-E_{l2}/T_{\text{ex}}}} \left[ \frac{1 - e^{-h\nu_1/k_B T_{\text{ex}}}}{1 - e^{-h\nu_2/k_B T_{\text{ex}}}} \right], \quad (5)$$

where $A_1$ and $A_2$ are the Einstein spontaneous emission coefficients, $g_{u1}$ and $g_{u2}$ are the statistical weights of the upper

levels, $E_{l1}(K)$ and $E_{l2}(K)$ are the excitation energies of the lower levels, and $Q_1$ ($N_1^{\text{tot}}$) and $Q_2$ ($N_2^{\text{tot}}$) are the partition functions (total column density) of the molecule/atom for line 1 and line 2, respectively. Therefore, if the two lines are from the same molecule/atom transitions, the opacity ratio is independent of the partition function and the total column density.

In the third case, the transition $\nu_1$ is optically thin, while the transition $\nu_2$ is optically thick. The integrated intensity ratio is

$$R_{\nu_1/\nu_2} = \frac{[J_{\nu_1}(T_{\text{ex}}) - J_{\nu_1}(T_{\text{bg}})]\tau_1}{[J_{\nu_2}(T_{\text{ex}}) - J_{\nu_2}(T_{\text{bg}})]\tau_2}. \quad (6)$$

In Figure 8, the left panels shows the results for the ratios CO(8-7)/CO(5-4) and CO(11-10)/CO(5-4) as a function of $T_{\text{ex}}$, assuming these two lines are both optically thin (solid line) and optically thick (dashed line). The latter case slowly changes or is likely constant for $T_{\text{ex}} \ge 40$ K. The horizontal lines show the observed ratios, which are CO(8-7)/CO(5-4) ~0.8 and CO(11-10)/CO(5-4) ~0.35. The range of $T_{\text{ex}}$ can be determined by the intersection of this line and these two extreme cases above (cyan-filled region), which is not far from the dust temperature in this region (≈30 K; see Roccatagliata et al. 2013). If we assume that these three CO lines trace the same excitation temperature, the $T_{\text{ex}}$ spans from 40 to 70 K. In Figure 8, the right panel shows the line ratio (contour lines) as a function of both $T_{\text{ex}}$ and $\tau$ by assuming CO(5-4) is optically thick and CO(11-10) and CO(8-7) are optically thin. For $T_{\text{ex}} \sim 40-70$ K, $\tau_{\text{CO(8-7)}} \sim 0.8$, and $\tau_{\text{CO(11-10)}} \sim 0.4$. We use $\tau_{\text{CO(11-10)}} \sim 0.4$ because it likely manifests our assumption to compute the total density of CO gas, which results in $N_{\text{CO}} \sim 2 \times 10^{16}-5 \times 10^{17} \text{ cm}^{-2}$. We rule out the case of CO(5-4) being optically thin, while CO(8-7) and CO(11-10) are optically thick (not shown here) because the optical depth returns to be greater than unity, which contradicts our assumption.

### 3.5.2. Non-LTE with RADEX

For the case of non-LTE conditions, we model the integrated intensity for each transition line of CO by the radiative transfer
<figure><img src="image_8.png" /><figcaption>**Figure 6.** PV plots of CO(5-4), CO(8-7), CO(11-10), and [O I] in contours enhanced with orange scale along the cuts shown in Figure 2. Five contours from 3 K to the peak value of the brightness temperature are plotted for each cut. The orange scale spans from 0 to 1.5 times the maximum temperature.</figcaption></figure>

program RADEX (van der Tak et al. 2007). To constrain $T_{\text{kin}}$, $n_{\text{H}_2}$, and $N_{\text{CO}}$, we vary these three parameters at the same time and fit the diagram of the integrated intensity (*I*) versus the quantum upper level (*J*). We use a Bayesian approach to

perform the fits. The parameter space is investigated using an ensemble of Markov chains following the affine-invariant Markov Chain Monte Carlo (MCMC; Goodman & Weare 2010) algorithm. In particular, inspired by Yang et al. (2017),
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Average spectra of four transition lines with the fitted Gaussians (red curves). The mean and standard deviation (of the main Gaussian) are noted in each panel. The yellow-filled area corresponds to the signal range, which lies within [$\mu - 3\sigma, \mu + 3\sigma$].</figcaption></figure>

we couple the emcee (Foreman-Mackey et al. 2013) library with a Python wrapper of RADEX called *pyradex*<fnref n="12" /> to obtain RADEX results for each iteration and sample the posterior probability distribution function (PDF) with emcee. We use a total of 400 walkers, iterated with 2000 “burn-in” steps followed by 10,000 main steps, to sample the posterior PDF.

The top panel of Figure 9 shows the best-fit model (orange line and yellow shaded area) to the observed data (black dots). The triangle plot in the bottom panel shows the distribution of the best-fit parameters. The vertical dashed lines mark the 16th and 84th percentiles, which indicate the lower and upper values of the uncertainties, respectively. The posterior PDF for the gas density is peaked at $n_{H_2} \simeq 1.5 \times 10^5 \text{ cm}^{-3}$ and $N(\text{CO}) \simeq 10^{16} \text{ cm}^{-2}$, while the kinetic temperature is not constrained, with $T_{\text{kin}}$ uncertainty ranging from $\approx 219$ to 724 K.

Our constraint on the CO column density of $10^{16}-5 \times 10^{17} \text{ cm}^{-2}$ and on the gas volume density of $10^5 \text{ cm}^{-3}$ are similar to the ones of the neighbor Treasure Chest globule ($N(\text{CO}) = 8.5 \times 10^{16}$

$-4.9 \times 10^{18} \text{ cm}^{-2}$ and $n_{H_2} \sim 10^5-10^6 \text{ cm}^{-3}$; see Mookerjea et al. 2019).

# 4. Discussion

## 4.1. *Motion of the G287.76-0.87 Pillar*

Determining the direction of motion relative to the illuminated source is important for determining the effects of stellar feedback on the pillar. To put the observed velocity of the [C II] emission into context, we first assess its relation to the stars that are the source of ionization and wind energy. The most luminous star in the region, $\eta$ Carina itself, has a velocity of $-19.7 \text{ km s}^{-1}$ with respect to LSR (Smith 2004). A survey of the radial velocities of 63 O-type stars in Carina has an average velocity of $-11.0 \text{ km s}^{-1}$, with a wide dispersion of $9 \text{ km s}^{-1}$ (Kiminki & Smith 2018). The pillar structure that is the subject of this paper is illuminated from the NW in the direction of Tr 16, which contains the most powerful Wolf-Rayet and O stars. The average velocity of O stars in Tr 16 is $-15.1 \text{ km s}^{-1}$ (Kiminki & Smith 2018). It is important to note the wide dispersion in velocities among stars in each cluster and differences between clusters of order $10 \text{ km s}^{-1}$.

<footnote n="12">¹² https://github.com/keflavich/pyradex</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Left: the ratio of the integrated intensity of CO(8-7)/CO(5-4) (top) and of CO(11-10)/CO(5-4) (bottom) as a function of the excitation temperature when they are both optically thin (solid line) and both optically thick (dashed line). Right: their integrated intensity ratios are computed as a function of the excitation temperature and the opacity, assuming that CO(5-4) is optically thick. Dashed-dotted lines mark the observed ratio calculated from the average integrated intensity of the lines.</figcaption></figure>

Now consider the molecular gas velocities. The MOPRA Southern Galactic Plane Survey shows a cold molecular gas traced by CO(1–0) emission at velocities around −20 km s⁻¹ (Rebolledo et al. 2016). The pillar G287.76-0.87 is a part of this complex of molecular gas that is currently a work surface exposed to stellar radiation and winds. The dispersion of the velocity of CO is ~6 km s⁻¹, so a core formed from this cold gas could have a velocity anywhere from −17 to −23 km s⁻¹.

If the pillar gas started at the same −20 km s⁻¹ as the cold molecular gas nearby and the pressure from Tr 16 arose at a slightly greater distance than the pillar, it would be pushed toward more negative velocities, $v < -20$ km s⁻¹. We measure a line-of-sight velocity of −14 km s⁻¹, which is contrary to the direction of expansion based on the assumptions of the previous sentence. However, the present-day cold molecular gas in the vicinity may not be representative of the initial conditions for the pillar because the portion of the cloud that comprises the pillar may have been part of the parent clump that formed Tr 16. If that is the case, then we should compare its velocity with that of the stars, and the acceleration would lead to $v < -15.1$ km s⁻¹.

Nevertheless, the observed radial velocity of the pillar still indicates that it is moving away from us by ~1 km s⁻¹ (averaged over the pillar). In addition, this is only the radial

component; the full 3D velocity depends on the angle between the ray from the pillar toward Tr 16 and the line of sight. This angle could be quite large because the illumination appears in the sky’s plane. Therefore, it is ambiguous to know whether the pillar is pushed away or toward us.

The Hubble image in Figure 10 shows a sharp, bright, exposed rim to the NW, with a dark core, qualitatively suggesting that the illumination source is approximately the same distance as the core but perhaps somewhat farther away so that a ray from Tr 16 to the molecular core lands on the NW side of the core just on the back side.

## 4.2. Stellar Feedback on the G287.76-0.87 Pillar

As discussed in Section 4.1, the interpretation of the G287.76-0.87 motion is uncertain. As for a demonstration, consider that the pillar is accelerated by 1 km s⁻¹; the kinetic energy required is $10^{45}$ erg due to a mass of order $10^2 M_\odot$. With a projected separation between the pillar and Tr 16 of 10 pc, and a pillar diameter of 0.7 pc, only 0.03% of the stellar output is intercepted by it. Thus, the energy output from the cluster would need to be of order $10^{48}$ erg over a spherical shell at 10 pc. The extant stars could produce such energies. Indeed, the pressure that a H II region (i.e., hot gas) powered by an illumination source is
<figure><img src="image_11.png" /><figcaption>**Figure 9.** RADEX modeling results with the MCMC sample of three free parameter spaces $n(\text{H}_2)$, $T_{\text{kin}}$, and $N(\text{CO})$. Top: orange curve shows the best RADEX and emcee fitting to the intensity of three CO lines, and yellow ones plot the 100 “good models” within the uncertainty range from the 16th to the 84th percentiles. Bottom: distribution of the best-fit parameters, in which the main diagonal shows the marginalized distributions of each parameter, with vertical dashed lines indicating the 16th and 84th percentiles, and the 2D plots show the correlation between each pair of parameters. The solid orange lines locate the best-fit model in the parameter space.</figcaption></figure>

$P_{\text{ext,HII}}/k_{\text{B}} = \sqrt{3Q/4\pi r^3\alpha_r} T_{\text{HII}}$ (see Equation (2) in Harper-Clark & Murray 2009) with $\alpha_r \simeq 4 \times 10^{-13}$ the recombination coefficient and $Q = 9 \times 10^{50} \text{s}^{-1}$ the ionizing photons per second. For a $T_{\text{HII}} = 10^4 \text{ K}$, $P_{\text{ext,HII}}/k_{\text{B}} \sim 10^6 \text{ K cm}^{-3}$. The equivalent energy of the cluster over a spherical shell at 10 pc is of the order $10^{49} \text{ erg}$.

The photons from the radiation source Tr 16 also compress the pillar with a pressure of $P_{\text{ext,rad}} = L_{\text{bol}}/(4\pi r^2 c)$, in the limit where all photons incident on the pillar are absorbed by it and the dust and gas are coupled (Tielens 1983). Adopting the bolometric luminosity of $\simeq 10^7 L_{\odot}$ (see, e.g., Smith 2006; Harper-Clark & Murray 2009), the radiation pressure exerted on the pillar is about the same as the hot gas pressure, $P_{\text{ext,rad}}/k_{\text{B}} \sim 10^6 \text{ K cm}^{-3}$.

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Hubble H$\alpha$ image of the pillar covering the same region as Figure 2. The main body of the pillar is dark, apparently in silhouette. The rim, and in particular the northwest edge (toward the luminous stars in Tr 16), is bright, indicating the pillar is closer to us than the illumination source.</figcaption></figure>

As the cloud contracts, the internal pressure due to the thermal gas acts outward, causing the cloud to reexpand. Taking into account the hydrostatic pressure by nonthermal gas motion, the internal gas pressure of $P_{\text{int,th}}/k_{\text{B}} = nT_{\text{eff}}$, with $T_{\text{eff}}$ as the effective gas temperature. As the velocity dispersion $\sigma = 0.92-1.18 \text{ km s}^{-1}$ (Figure 7), the representative effective temperature $T_{\text{eff}}^{\text{average}} = (0.92^2 - 1.18^2) \times 10^{10} (\mu \text{mH}/k_{\text{B}})$ K with $\mu = 2.8$; thus, the average internal gas pressure is $P_{\text{int,th}}^{\text{average}}/k_{\text{B}} \sim 10^7 \text{ K cm}^{-3}$. As a consistency check, the thermal pressure of the gas could be estimated via the empirical relation to the radiation field as $P_{\text{int,th}}^{\text{PDR}}/k_{\text{B}} = 2.1 \times 10^4 G_{\text{UV}}^{0.9}$ derived from the PDR models (Wu et al. 2018). For $G_{\text{UV}} \sim 1000 \times G_0$ at the location of our pillar (Roccatagliata et al. 2013), $P_{\text{int,th}}^{\text{PDR,average}}/k_{\text{B}} \sim 10^7 \text{ K cm}^{-3}$. This internal pressure is about an order of magnitude higher than the external pressure of hot gas and radiation. This could be due to (1) the overestimation of $T_{\text{eff}}$ (an average over integrated spectrum that could not represent the local variation) and (2) the gas heated externally and internally. Nevertheless, this simply implies that the pillar is supported against the feedback from these external pressures, which seems consistent with a small relative velocity with respect to $\eta$ Carina and is plausible for the strong turbulent motion in the tail as seen in [O I] as the gas is leaked through it (the gas density is lower at the tail).

Interestingly, the stellar wind pressure at 10 pc is $P_w/k_B \simeq 10^9 \text{ K cm}^{-3}$ at 10 pc assuming a standard mass-loss rate of $\sim 10^{-3} M_{\odot} \text{ yr}^{-1}$ (see Harper-Clark & Murray 2009), which is much higher than all above pressures. This implies that the stellar wind either escaped (leaked) through the holes (lowest density) to sculpt the more distant interstellar medium, or was dissipated, or the wind luminosity is overestimated. In reality, a combination of these three factors might be the case. The former is phenomenally similar to the expansion of gas around the Orion Nebula Cluster (Pabst et al. 2019) or 30 Doradus (Melnick et al. 2021) and could be plausible for the
giant bubble around Tr 16. The latter phenomenon could be
due to the shocked gas in which the stellar wind loses its kinetic
energy. A more systematic study of the velocities of the Carina
pillars could yield a real determination of the feedback of the
star cluster onto massive cores near one of the most powerful
stellar clusters of the Milky Way.

5. Summary

In this paper, we studied the kinematics of the pillar G287.76-0.87 in the southern region of the Carina Nebula, using the observational data from SOFIA/GREAT, including CO(5-4), CO(8-7), CO(11-10), and [O I] 63 µm. Our main findings are summarized as follows:

1. We analyzed the channel maps and the resolved spectral lines and showed that the line of sight velocity of the pillar is −14 km s<sup>−1</sup>, which is more redshifted than the components of the Tr 16 cluster. The optical image reveals that the rim and the edge toward the luminous stars in Tr 16 are bright, indicating that the pillar is closer to us than the illumination source. Thus, the LOS velocity contradiction might be solely due to the project effect.

2. We made the moment maps of the four lines. The zeroth moment maps show that the pillar's northwest region (toward $\eta$ Carina) shows strong emission because of the radiative excitation and/or interaction with the H II expansion shell. The [O I] 63 $\mu$m line partially correlates very well with the PAH continuum emission, and CO lines rise behind the atomic line, which illustrates the structure of the PDR. Second-moment maps show that the motion in the tail is more turbulent.

3. We extracted the PV diagram along several directions and indicated that the H II expansion shell influences the gas motion in the northwest-to-southeast direction (Cut 1 in Figure 2), while the gas is less dynamic in the perpendicular direction (cut 2 in Figure 2). The diagram through the southernmost (cut 3 in Figure 2) indicates a diffuse medium.

4. We constrained the physical properties of the G287.76-0.87 pillar under LTE and non-LTE conditions. For LTE, our method estimated the CO column density of $N_{CO} \sim 2 \times 10^{16} - 5 \times 10^{17}$ cm$^{-2}$. For non-LTE, we used RADEX with MCMC to estimate the gas number density, kinetic temperature, and CO column density. These estimations yield $n_{H_2} \simeq 1.5 \times 10^5$ cm$^{-3}$ and $N_{CO} \simeq 10^{16}$ cm$^{-2}$.

5. We discussed that the internal pressure of the gas within the G287.76-0.87 is sufficiently high that this pillar is hard to be pushed by the external hot gas and radiation pressures from Tr 16.

Acknowledgments

Based on observations made with the NASA/DLR Strato-
spheric Observatory for Infrared Astronomy (SOFIA). SOFIA is
jointly operated by the Universities Space Research Association,
Inc. (USRA), under NASA contract NNA17BF53C, and the
Deutsches SOFIA Institut (DSI) under DLR contract 50 OK
0901 to the University of Stuttgart. N.D.T acknowledges
financial support from the French National Research Agency
(ANR) through the Programme d'Investissement d'Avenir under
contracts ANR-11-LABX-0058_NIE and ANR-17-EURE-0024
within the Programme d'Investissement d'Avenir ANR-10-
IDEX-0002-02, and from the European Research Council
(ERC) via the ERC Synergy Grant ECOGAL (grant 855130).

*Software:* GILDAS/CLASS (Pety 2005; GILDAS Team 2013), APLpy (Robitaille & Bressert 2012), RADEX (van der Tak et al. 2007), pvextractor (Ginsburg et al. 2016), kalibrate (Guan et al. 2012), emcee (Foreman-Mackey et al. 2013), pyradex (https://github.com/keflavich/pyradex).

Appendix A
Channel Maps of CO

Figures 11–13 show the channel maps of the ¹²CO(5-4), ¹²CO(8-7), and ¹²CO(11-10). Owing to the larger beam sizes, ¹²CO maps are able to spatially resolve the pillar structure. However, the bulk velocity is also ≃ −14 km s⁻¹ as seen in [O I] (see Figure 3).
<figure><img src="image_13.png" /><figcaption>Figure 11. Channel maps of CO(5-4).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 12. Channel maps of CO(8-7).</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure 13. Channel maps of CO(10-11).</figcaption></figure>

## Appendix B
### Moment Maps of CO Transition Lines

The left panels in Figure 14 show the moment-zero maps of $^{12}$CO(5-4), $^{12}$CO(8-7), and $^{12}$CO(11-10), respectively, from top

to bottom. Similarly, the right panels show the moment-one maps. These $^{12}$CO lines are bright in emission, but their maps are unable to resolve the pillar structure.
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Integrated intensity maps (left panel) and mean-velocity maps (right panel) of CO(5-4), CO(8-7), and CO(11-10) lines. The contours show the level of the integrated intensity.</figcaption></figure>

ORCID iDs

Ngo-Duy Tung https://orcid.org/0009-0009-8545-2682

Le Ngoc Tram https://orcid.org/0000-0002-6488-8227

Archana Soam https://orcid.org/0000-0002-6386-2906

William T. Reach https://orcid.org/0000-0001-8362-4094

Ed Chambers https://orcid.org/0000-0003-4195-1032
Blesson Mathew https://orcid.org/0000-0002-7254-191X

Heiko Richter https://orcid.org/0000-0002-3639-6997

References

Anderson, L. D., Makai, Z., Luisi, M., et al. 2019, ApJ, 882, 11

Durán, C. A., Güsten, R., Risacher, C., et al. 2021, ITTST, 11, 194

Elmegreen, B. G., & Lada, C. J. 1977, ApJ, 214, 725

Foreman-Mackey, D., Hogg, D. W., Lang, D., et al. 2013, PASP, 125, 306

Ginsburg, A., Robitaille, T., & Beaumont, C. 2016, pvextractor: Position-Velocity Diagram Extractor, Astrophysics Source Code Library, ascl:1608.010

Goldsmith, P. F. 2019, ApJ, 887, 54

Goodman, J., & Weare, J. 2010, CAMCS, 5, 65

GILDAS Team 2013, GILDAS: Grenoble Image and Line Data Analysis Software, Astrophysics Source Code Library, ascl:1305.010

Guan, X., Stutzki, J., Graf, U. U., et al. 2012, A&A, 542, L4

Harper-Clark, E., & Murray, N. 2009, ApJ, 693, 1696

Hester, J. J., Scowen, P. A., Sankrit, R., et al. 1996, AJ, 111, 2349

Hollenbach, D. J., & Tielens, A. G. G. M. 1997, ARA&A, 35, 179

Hsieh, B.-C., Wang, W.-H., Hsieh, C.-C., et al. 2012, ApJS, 203, 23

Kahn, F. D. 1954, BAN, 12, 187

Kavak, Ü., Bally, J., Goicoechea, J. R., et al. 2022, A&A, 663, A117

Kiminki, M. M., & Smith, N. 2018, MNRAS, 477, 2068

Kirsanova, M. S., Ossenkopf-Okada, V., Anderson, L. D., et al. 2020, MNRAS, 497, 2651

Kirsanova, M. S., Pavlyuchenkov, Y. N., Olofsson, A. O. H., Semenov, D. A., & Punanova, A. F. 2023, MNRAS, 520, 751

Melnick, J., Tenorio-Tagle, G., & Telles, E. 2021, A&A, 649, A175

Mookerjea, B., Sandell, G., Güsten, R., et al. 2019, A&A, 626, A131

Pabst, C., Higgins, R., Goicoechea, J. R., et al. 2019, Natur, 565, 618

Pety, J. 2005, in SF2A-2005: Semaine de l'Astrophysique Francaise, ed. F. Casoli et al. (Les Ulis: EdP-Sciences), 721

Rathborne, J. M., Brooks, K. J., Burton, M. G., Cohen, M., & Bontemps, S. 2004, A&A, 418, 563

Rebolledo, D., Burton, M., Green, A., et al. 2016, MNRAS, 456, 2406

Risacher, C., Güsten, R., Stutzki, J., et al. 2018, JAI, 7, 1840014

Robitaille, T., & Bressert, E. 2012, APLpy: Astronomical Plotting Library in Python, Astrophysics Source Code Library, ascl:1208.017

Roccatagliata, V., Preibisch, T., Ratzka, T., & Gaczkowski, B. 2013, A&A, 554, A6

Schneider, N., Güsten, R., Tremblin, P., et al. 2012, A&A, 542, L18

Smith, N. 2006, MNRAS, 367, 763

Smith, N. 2004, MNRAS, 351, L15

Smith, N., Stassun, K. G., & Bally, J. 2005, AJ, 129, 888

Strömgren, B. 1939, ApJ, 89, 526

Thompson, M. A., White, G. J., Morgan, L. K., et al. 2004, A&A, 414, 1017

Tielens, A. G. G. M. 1983, A&A, 119, 177

van der Tak, F. F. S., Black, J. H., Schöier, F. L., et al. 2007, A&A, 468, 627

White, G. J. 1993, A&A, 274, L33

White, G. J., & Gee, G. 1986, A&A, 156, 301

Wu, R., Bron, E., Onaka, T., et al. 2018, A&A, 618, A53

Yang, C., Omont, A., Beelen, A., et al. 2017, A&A, 608, A144

Young, E. T., Becklin, E. E., Marcum, P. M., et al. 2012, ApJL, 749, 17