Research Article NAO, North Atlantic Oscillation; RPC, ratio of predictable
components. Daniel J. Brener, Higgs Centre for Theoretical Physics, School of
Physics & Astronomy, The University of Edinburgh, James Clerk Maxwell Building
Room 4301, Kings Buildings, Peter Guthrie Tait Road, Edinburgh, EH9 3FD, UK
<EMAIL_ADDRESS>Science and Technology Facilities Council, UK.

# A possible link between ergodicity and the signal-to-noise paradox

Daniel J. Brener

###### Abstract

This short letter raises the possibility that ergodicity concerns might have
some bearing on the signal-to-noise paradox. This is explored by simply
applying the ergodic theorem of statistical mechanics to the theory behind
ensemble weather forecasting to find a source of error in the ensemble mean.
The problem of non-ergodic systems and models in the case of weather
forecasting is discussed.

###### keywords:

atmosphere, ergodicity, ensembles, NAO, predictability, signal-to-noise
paradox, climate

##  1 Introduction

The equations governing the atmosphere were shown in the seminal papers of
[12] and [13] as exhibiting sensitive dependence on their initial conditions,
commonly known as chaos. Since it is not possible to know the exact initial
conditions for any event, there will exist a limit on the atmosphere’s
predictability, even with a theoretically perfect model. In order to account
for this inherent uncertainty the method of ensemble weather prediction was
developed [6]. This process uses multiple realisations of the same numerical
model, with each realisation differing by tiny perturbations in the initial
conditions. Due to the chaotic nature of the equations these small
perturbations give rise to divergent predictions. The average value of an
observable, such as temperature, over several ensemble members, is interpreted
to be the best estimate of the said state of the observable at or over a
particular time. The standard deviation over the ensemble hence gives an
indicator of the uncertainty and chaotic state of the atmosphere, despite the
underlying equations governing the models being deterministic [14] [17].

The [22] review demonstrated a paradoxical result: that atmosphere-ocean
coupled climate and long range prediction ensemble models are better at
predicting reality, than they are at predicting themselves. They examined a
statistic derived from a Baysian framework for comparing ensemble models to
observations [23] known as the “ratio of predictable components" (RPC) for
several atmospheric parameters. The RPC can be written,

$RPC^{2}=\frac{r_{\bar{E}O}^{2}}{r_{\bar{E}E_{i}}^{2}},$ (1)

where $r_{\bar{E}O}$ is the correlation between the model ensemble mean and
the observations, and $r_{\bar{E}E_{i}}$ is the correlation of the ensemble
mean and any single ensemble member. It was found that the correlation between
the ensemble mean and observations is often much greater than the correlation
between the ensemble mean and any single ensemble member ($RPC>1$). Hence, a
signal-to-noise paradox - the model predicts reality better than it predicts
itself. This was first raised as possibility by [10] and has since been
reported by many different groups such as [5], [20] [24], [3], [29] and most
recently [4].

Over the past few years a great deal of effort from the community has been put
into finding different resolutions to the paradox, and different scenarios
where it arises. [26] showed that for the North Atlantic Oscillation (NAO),
the paradox may be interpreted as simply being a probabilistically under
confident forecast for occurrences of high magnitude NAO. A possible
explanation for this is that the model has a reduced persistence of particular
atmospheric regimes, especially in the Northern Hemisphere [27] [25]. In terms
of model dynamics, this might be due to weak atmosphere eddy feedback [21],
[7], or weak ocean-atmosphere coupling in models [16]. The early ideas around
the paradox were that through further investigation, potentially direct
improvements to the physics of the models could be found which would increase
the correlation between the model and the ensemble mean. This would then, in
theory, increase the overall forecast skill. So far, such necessary
improvements have yet to be discovered and implemented. Taking a slightly
different angle, [2] argue that in some cases, the signal-to-noise paradox
should not in fact be considered paradoxical due to the assumption that the
forecast error is related to the correlation of the ensemble mean with the
observations, which is not necessarily always the case. It is in this more
statistical direction that this Letter takes us.

In this Letter we will explore the idea that part of the paradox might be due
to issues with the accuracy of the ensemble mean with respect to representing
the true average. In statistical mechanics there is the fundamental principle
that macrostates emerge as ensemble averages of microstates. This is
essentially the idea being exploited by ensemble numerical weather prediction.
By running multiple models with different initial conditions one plays out
alternative realities in the systems phase space, the average over which is
taken to be the most likely macrostate the atmosphere will evolve to. The
system phase space for a theoretically perfect model should be identical to
the real atmosphere’s phase space. However, we must rely on theoretically
imperfect models, therefore the phase space explored is some combination of
the phase space inherent to the model and that of the atmosphere as given by
the observations.

If we are to consider each ensemble member, and the observations as equally
likely microstates in some probability space (e.g. a NAO probability
distribution function) then it is necessary that the Ergodic Theorem hold for
statistical averaging. If this theorem does not always hold then direct
comparison between the model ensemble mean and single members, or observations
is not valid, as the mean has ceased to be representative of the phase space.
The purpose of this short Letter is to theoretically illustrate this
ergodicity assumption. The ideas presented here are not intended to resolve
the paradox entirely as there are likely other contributing issues related to
both modelling and the definitions used in the computation of the RPC. (e.g.
[2], [30], [8] and [15]) We hope this paper will act as a catalyst for others
to consider ways in which the ergodicity assumption might be tested more
rigorously, and as a reminder of it.

##  2 Ergodic theory

###  2.1 The ergodic theorem

We will informally introduce the theorem but for a formal approach see [9].
For some observable $f(x)$, whose evolution in time is governed by its phase-
space distribution function $P(x,t)$, the average over an interval $\Delta
T=T-T_{0}$ can be found by,

$\bar{f}(x)=\lim_{\Delta T\rightarrow\infty}\frac{1}{\Delta
T}\int_{T_{0}}^{T}f(x)P(x,t)dt=\int\mu f(x)dx,$ (2)

where $\mu$ is known as the invariant measure for the phase-space $x$. This is
also known as the ergodic hypothesis and was proved by [1]. In words, it tells
us that the long time average of an observable is equivalent to the phase
space average. The theorem is valid provided the existence of both the
invariant measure $\mu$ and that $P$ is stationary over the averaging
interval. Ergodicity is often a necessary assumption when modelling physical
systems. However an exact equivalence between time average and ensemble mean
is only possible when the system is allowed to visit all the possible
microstates over a long period of time. For non-equillibirum systems this, in
general, cannot be guaranteed.

###  2.2 Application to ensemble modelling

Next we will introduce how the ensemble mean is generally defined in ensemble
weather models, and then show how this requires an assumption of ergodicity.
The ensemble average, $E(f(t))$ for a observable $f$ at a time $t=t_{0}+\delta
t$ where $t_{0}$ is the model start time is,

$\begin{split}E(f(t))=\frac{1}{n}\sum_{i=1}^{n}f(t)_{i}\end{split}$ (3)

where $f(t)_{i}$ is the observable from the $i^{th}$ member of the ensemble of
size $n$. In the language of statistical mechanics, we have taken an average
over all the microstates (ensemble members) in the system phase space at an
instant $t$, who’s size is defined by $n$. The significance of the ensemble
size, $n$ has been well studied on all temporal and spatial scales [17], [11].
The forecast skill grows with $n$ due to the suppression of unpredictable
noise.

Now consider each ensemble member as representing a microstate of the model
phase space, with distribution function $P_{i}(t)$. Then the average for some
observable $f$ over an interval $\Delta T$ is,

$\begin{split}E(f(\Delta T))&=\frac{1}{\Delta
T}\int_{T_{0}}^{T}{f_{1}(t)P_{1}(t)+f_{2}(t)P_{2}(t)+...+f_{n}(t)P_{n}(t)dt}\end{split}$
(4)

where the ensemble, or number of microstates, is size $n$. Mathematically, the
real phase space is not discrete, so ideally we need a ensemble of infinite
size. If the ergodic theorem is to hold, then the average must be taken over a
long period of time. Applying both these constraints to equation 4 give us the
“real ensemble average", $E_{R}$

$E_{R}=\lim_{\begin{subarray}{c}\Delta T\rightarrow\infty\\\
n\rightarrow\infty\end{subarray}}\left(E(f(\Delta T))\right),$ (5)

so we can see there are several assumptions which go into the ensemble mean
arising from ergodicity, the most significant of which is that the
distribution functions for each ensemble member are identical and stationary.
If $P(t)$ is not stationary, and allowed to vary over the averaging period
$\Delta T$, an error is introduced into the ensemble mean which is not
necessarily accounted for by the spread of the ensemble members.

The first obvious mitigation strategy is to simply increase the ensemble size.
However this has been shown by [22] to have no relation to the signal-to-noise
paradox. There is weak evidence that the paradox exists on longer than multi-
decadal timescales [22]. The reason for this may be understood from an ergodic
perspective, as to define a mean climate state, the interval $T-T_{0}$ over
which we compute the average, is sufficiently long that the distributions
$P_{i}$ converge [28]. A second idea might be to increase the averaging window
$\Delta T$, it would be interesting to see what effect this has on the RPC,
however one cannot increase it too much else the definition of the observable
would be changed. In other words, for very large averaging windows, the
observable simply becomes a different kind of climatic average: sub-seasonal
to seasonal, seasonal to annual etc.

Our key point is that in a pure sense one cannot assume that the underlying
probability distribution function is the same for each member, even though
they are constructed from the same model. Each ensemble member represents,
after a short time, a unique “reality" as each of them eventually forgets
their initial conditions. To improve the ensemble spread and hence the
forecast skill further, stochastic physics is used. However, this has
implications for computing and interpreting the ensemble mean. One cannot
factorise a single $P_{i}(t)$ in the integrand of equation 4. This is an
assumption underlying the ensemble mean - that all the distributions are
identical. This should be an acceptable approximation to make for many
atmospheric processes, but not all. For seasonal predictions it is possible
that the models are more sensitive to this assumption due to the exponential
error growth, and hence faster divergence of the probability density
functions. Stating this in an alternative, more practical way: the Ergodic
Theorem allows us to make the approximation in equation 5 if the observable
$f$ follows a stationary distribution. This means that the underlying
probability distribution $P(t)$, of the observable must not change between
different ensemble members at each instant over the averaging interval
$T-T_{0}$. In the limit $t\rightarrow\infty$, each member will forget its
initial conditions, but for long range forecasts, one must preserve the key
invariants for the observable: its probability distribution function must
remain stationary with respect to the averaging window and the invariant
measure must be invariant.

###  2.3 Discussion

One way to investigate this idea would be to develop an “ensemble of
ensembles". We borrow this idea from a technique for exploring a phase space
more completely in direct numerical simulation of the Navier-Stokes equations.
At each timestep a small perturbation is made to each member of the ensemble
from which several different realisations are computed - the difference
between them generating a finite time Lyapunov exponent in the usual way. One
could investigate ergodicity by at each time step generating a new ensemble
with enough members to compute distributions, such as the Lyapunov spectrum.
If the spectrum varies between the ensemble members then this would imply that
in a fundamental way, the different ensemble members evolve into different
state spaces. Thus interpreting the ensemble mean to be the most likely
weather scenario would not be justified.

One approach to relax the ergodicity constraint but still make use of
ensembles is to enforce dynamical invariants. These are quantities that are
conserved along the phase space path of the system. For example, in the
probability theory for non-equilibrium gravitational systems developed by
[18], dynamical invariants are employed with the explicit purpose of not
relying ergodicity assumptions. For turbulent fluid systems, the most obvious
dynamical invariant would be the Lyapunov spectrum and fractal dimension,
besides the more common helicity quantities, such as vorticity.

Recently [19] used the same principle from ergodic theory to developed a
weather forecasting machine learning training method. Their algorithm enforced
the dynamical invariant of the Lyapunov spectrum and fractal dimension. They
tested their neural network on the classic Lorenz 1996 chaotic dynamical
system and the two-layer baroclinic model of Charney and Strauss with added
baroclinic instabilities. Including the ergodic constraints improved the
forecast skill of the neural network. We speculate that the RPC value will be
improved for long-range forecasts if these invariant constraints can be
carefully applied to the ensembles. Each observable should have its own
invariant measure for the ergodic theorem to hold for the ensemble. However,
applying such an invariant might also constrain the ensemble in an unhelpful
way, as it is possible that the real world observable does not obey the
invariant anyway.

The link to the paradox is not that there is something special about the
observations versus the model, but that the definition of the ensemble mean
does not represent the true microstate of the ensemble. This connects directly
to the idea of [2], that the error in the ensemble mean results in a
meaningless RPC. As in general, the correlation between ensemble mean and
ensemble members does not equal the correlation between ensemble mean and
observations. Our Letter therefore is pointing to a possible source of the
error in the ensemble mean.

##  3 Concluding remarks

The core idea in this Letter is that imperfectly constrained ergodicity will
introduce a source of error to the ensemble mean, and therefore the signal-to-
noise paradox may naturally arise in certain cases. If we are to treat every
ensemble member as an equally likely weather scenario or in statistical
physics terms, a microstate, then each ensemble member can evolve according to
its own probability distribution function which must remain stationary over
the averaging window and stationary relative to all the other ensemble
members. Effectively each ensemble member is capable of introducing its own
inherent source of error to the mean. Likewise ergodicity may be broken if the
ensemble members or the observable in question, do not have or obey an
invariant measure. We do not suggest that these divergences from ergodicity
are necessarily a common occurrence, or that they entirely resolve the
paradox. There are already many different suggested reasons why a signal-to-
noise paradox can arise, the purpose of this Letter is simply to highlight
another consideration. We highlight two potential ways of investigating the
role of ergodicity within the paradox: examining the finite time Lyapunov
exponent evolution for a set of ensemble members and constraining the
dynamical invariants such as the Lyapunov spectrum for a complete ensemble.

## ACKNOWLEDGEMENTS

D.J.B. is funded by the UK Science and Technology Facilities Council. He is
also particularly grateful to Professor Adam Scaife who introduced him to the
problem, Dr Jose Rodriguez, and Dr Leon Hermanson of the Met Office and
Professor Jorge Peñarrubia of The University of Edinburgh, for giving their
time to insightful tête-à-tête discussions. He also thanks Alexander Tully of
The University of Liverpool for comments on the manuscript. For the purpose of
open access, the author has applied a Creative Commons Attribution (CC BY)
licence to any Author Accepted Manuscript version arising from this
submission.

## CONFLICT OF INTEREST

The author declares no conflict of interest.

## References

  * Birkhoff [1931] Birkhoff, G. D. (1931) Proof of the ergodic theorem.  Proceedings of the National Academy of Sciences of the United States of America, 17, 656–660.  URL: http://www.jstor.org/stable/86016.
  * Bröcker et al. [2023] Bröcker, J., Charlton–Perez, A. J. and Weisheimer, A. (2023) A statistical perspective on the signal-to-noise paradox.  Quarterly Journal of the Royal Meteorological Society, 149, 911–923.  URL: https://rmets.onlinelibrary.wiley.com/doi/abs/10.1002/qj.4440.
  * Charlton-Perez et al. [2019] Charlton-Perez, A., Bröcker, J., Stockdale, T. N. and Johnson, S. J. (2019) When and where do ecmwf seasonal forecast systems exhibit anomalously low signal-to-noise ratio?  Quarterly Journal of the Royal Meteorological Society, 145, 3466–3478.  URL: https://doi.org/10.1002/qj.3631.
  * Dunstone et al. [2023] Dunstone, N., Smith, D. M., Hardiman, S. C., Hermanson, L., Ineson, S., Kay, G., Li, C., Lockwood, J. F., Scaife, A. A., Thornton, H., Ting, M. and Wang, L. (2023) Skilful predictions of the summer north atlantic oscillation.  Communications Earth & Environment, 4, 409.  URL: https://doi.org/10.1038/s43247-023-01063-2.
  * Eade et al. [2014] Eade, R., Smith, D., Scaife, A. A., Wallace, E., Dunstone, N., Hermanson, L. and Robinson, N. (2014) Do seasonal-to-decadal climate predictions underestimate the predictability of the real world?  Geophysical research letters, 41, 5620–5628.  URL: https://doi.org/10.1002/2014gl061146.
  * Epstein [1969] Epstein, E. S. (1969) Stochastic dynamic prediction.  Tellus, 21, 739–759.
  * Hardiman et al. [2022] Hardiman, S. C., Dunstone, N. J., Scaife, A. A., Smith, D. M., Comer, R., Nie, Y. and Ren, H.-L. (2022) Missing eddy feedback may explain weak signal-to-noise ratios in climate predictions.  npj Climate and Atmospheric Science, 5, 57.
  * Knight et al. [2022] Knight, J. R., Scaife, A. A. and Maidens, A. (2022) An extratropical contribution to the signal-to-noise paradox in seasonal climate prediction.  Geophysical Research Letters, 49, e2022GL100471.  URL: https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2022GL100471.  E2022GL100471 2022GL100471.
  * Kong [2019] Kong, T. (2019) Ergodic theory, entropy and application to statistical mechanics.  URL: https://api.semanticscholar.org/CorpusID:208176466.
  * Kumar [2009] Kumar, A. (2009) Finite samples and uncertainty estimates for skill measures for seasonal prediction.  Monthly Weather Review, 137, 2622–2631.
  * Leutbecher [2019] Leutbecher, M. (2019) Ensemble size: How suboptimal is less than infinity?  Quarterly Journal of the Royal Meteorological Society, 145, 107–128.  URL: https://rmets.onlinelibrary.wiley.com/doi/abs/10.1002/qj.3387.
  * Lorenz [1963] Lorenz, E. N. (1963) Deterministic nonperiodic flow.  Journal of atmospheric sciences, 20, 130–141.
  * Lorenz [1969] — (1969) The predictability of a flow which possesses many scales of motion.  Tellus, 21, 289–307.
  * Murphy and Palmer [1986] Murphy, J. and Palmer, T. (1986) Experimental monthly long-range forecasts for the united-kingdom. 2. a real-time long-range forecast by an ensemble of numerical integrations.  Meteorological Magazine, 115.
  * O’Reilly et al. [2019] O’Reilly, C. H., Weisheimer, A., Woollings, T., Gray, L. J. and MacLeod, D. (2019) The importance of stratospheric initial conditions for winter north atlantic oscillation predictability and implications for the signal-to-noise paradox.  Quarterly Journal of the Royal Meteorological Society, 145, 131–146.  URL: https://rmets.onlinelibrary.wiley.com/doi/abs/10.1002/qj.3413.
  * Ossó et al. [2020] Ossó, A., Sutton, R., Shaffrey, L. and Dong, B. (2020) Development, amplification, and decay of atlantic/european summer weather patterns linked to spring north atlantic sea surface temperatures.  Journal of Climate, 33, 5939–5951.
  * Palmer and Hagedorn [2006] Palmer, T. and Hagedorn, R. (2006) Predictability of Weather and Climate.  Cambridge University Press.
  * Peñarrubia [2015] Peñarrubia, J. (2015) A probability theory for non-equilibrium gravitational systems.  Monthly Notices of the Royal Astronomical Society, 451, 3537–3550.  URL: https://doi.org/10.1093/mnras/stv1146.
  * Platt et al. [2023] Platt, J. A., Penny, S. G., Smith, T. A., Chen, T.-C. and Abarbanel, H. D. I. (2023) Constraining chaos: Enforcing dynamical invariants in the training of reservoir computers.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 33, 103107.  URL: https://doi.org/10.1063/5.0156999.
  * Scaife et al. [2014] Scaife, A. A., Arribas, A., Blockley, E. W., Brookshaw, A., Clark, R. T., Dunstone, N., Eade, R., Fereday, D., Folland, C. K., Gordon, M., Hermanson, L., Knight, J., Lea, D. J., MacLachlan, C., Maidens, A., Martin, M., Peterson, A., Smith, D., Vellinga, M., Wallace, E., Waters, J. and Williams, A. J. (2014) Skillful long-range prediction of european and north american winters.  Geophysical Research Letters, 41, 2514–2519.  URL: https://doi.org/10.1002/2014gl059637.
  * Scaife et al. [2019] Scaife, A. A., Camp, J., Comer, R., Davis, P., Dunstone, N., Gordon, M., MacLachlan, C., Martin, N., Nie, Y., Ren, H.-L. et al. (2019) Does increased atmospheric resolution improve seasonal climate predictions?  Atmospheric Science Letters, 20, e922.
  * Scaife and Smith [2018] Scaife, A. A. and Smith, D. (2018) A signal-to-noise paradox in climate science.  npj Climate and Atmospheric Science, 1, 28.  URL: https://doi.org/10.1038/s41612-018-0038-4.
  * Siegert et al. [2015] Siegert, S., Stephenson, D. B., Sansom, P. G., Scaife, A. A., Eade, R. and Arribas, A. (2015) A Bayesian Framework for Verification and Recalibration of Ensemble Forecasts: How Uncertain is NAO Predictability?  Journal of Climate, 29, 995–1012.  URL: https://doi.org/10.1175/JCLI-D-15-0196.1.
  * Stockdale et al. [2015] Stockdale, T. N., Molteni, F. and Ferranti, L. (2015) Atmospheric initial conditions and the predictability of the arctic oscillation.  Geophysical Research Letters, 42, 1173–1179.  URL: https://doi.org/10.1002/2014gl062681.
  * Strommen [2020] Strommen, K. (2020) Jet latitude regimes and the predictability of the north atlantic oscillation.  Quarterly Journal of the Royal Meteorological Society, 146, 2368–2391.
  * Strommen et al. [2023] Strommen, K., MacRae, M. and Christensen, H. (2023) On the relationship between reliability diagrams and the “signal-to-noise paradox”.  Geophysical Research Letters, 50, e2023GL103710.  URL: https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2023GL103710.  E2023GL103710 2023GL103710.
  * Strommen and Palmer [2018] Strommen, K. and Palmer, T. (2018) Signal and noise in regime systems: A hypothesis on the predictability of the north atlantic oscillation.  Quarterly Journal of the Royal Meteorological Society, 145, 147–163.  URL: https://doi.org/10.1002/qj.3414.
  * Tantet [2016] Tantet, A. J. J. (2016) Ergodic theory of climate: variability, stability and response.  Ph.D. thesis, Utrecht University.
  * Weisheimer et al. [2019] Weisheimer, A., Decremer, D., MacLeod, D., O’Reilly, C. H., Stockdale, T. N., Johnson, S. J. and Palmer, T. (2019) How confident are predictability estimates of the winter north atlantic oscillation.  Quarterly Journal of the Royal Meteorological Society, 145, 140–159.  URL: https://doi.org/10.1002/qj.3446.
  * Zhang [2019] Zhang, Wei; Kirtman, B. P. (2019) Understanding the signal-to-noise paradox with a simple markov model.  Geophysical Research Letters, 46, 13308–13317.  URL: https://doi.org/10.1029/2019gl085159.


# Star Formation Triggered by the expanding bubble S111

Bhaswati Mookerjea Tata Institute of Fundamental Research,
Homi Bhabha Road, Mumbai 400005, India

###### Abstract

This paper investigates the impact of radiative and mechanical feedback from
O-type stars on their parent molecular clouds and the triggering of formation
of future generation of stars. We study the infrared bubble S111 created by
the embedded massive stellar cluster G316.80-0.05. A significant fraction of
gas in shells created due to the compression of the ambient medium by
expanding bubbles is photodissociated by the stellar radiation. The kinematics
of the shells are thus best studied using spectroscopic observations of singly
ionized carbon, the most dominant species. We have used the velocity-resolved
maps of the ${}^{2}{\rm P}_{3/2}\rightarrow^{2}{\rm P}_{1/2}$ transition of [C
ii] at 158 µm, the $J$=2–1 transition of 13CO and C18O, and the $J$=1–0
transition of HCO+ to study the rim of the bubble S111 that partly coincides
with the southern part of the infrared dark ridge G316.75. The [C ii] spectra
conclusively show evidence of a shell expanding with a moderate velocity of
$\sim 7$ km s-1, which amounts to a kinetic energy that is $\sim$0.5–40 times
the thermal energy of the H ii region. The pressure causing the expansion of
the H ii region arises mainly from the hydrogen ionization and the dust-
processed radiation. Among the far-infrared sources located in the compressed
shells, we find the core G316.7799-0.0942 to show broad spectral features
consistent with outflow activity and conclude that it is a site of active star
formation. Based on the age of the H ii region we conclude that this expanding
H ii region is responsible for the triggering of the current star formation
activity in the region.

ISM: Bubbles – infrared: ISM – ISM: lines and bands –(ISM:) photon-dominated
region (PDR) –ISM: individual (S111) –ISM: individual (G316.75) – ISM:
kinematics and dynamics

††journal: APJ

##  1 Introduction

Massive stars play a dominant role in the injection of UV radiation into the
interstellar medium (ISM) and of mechanical energy through stellar winds and
finally through supernova explosions. The photoionization, radiation pressure,
and stellar winds from young O-type stars pump both energy and momentum into
the ambient material and disrupt the parent molecular cloud by creating H ii
regions and parsec-size bubbles enclosed by shells of denser swept-up material
(Churchwell et al., 2006, and references therein). It is envisaged that such
radiative and mechanical feedback from the massive stars regulate the star
forming environment locally as well as influence the evolution of the ISM in
the galaxies as a whole (Krumholz et al., 2018). Inventories of such bubbles
have been primarily derived from large-scale mid-infrared continuum images
observed using the Spitzer and AKARI telescopes (Deharveng et al., 2010;
Hanaoka et al., 2020).

The triple bubble system consisting of S109, S110, and S111 (Churchwell et
al., 2006) is created by the bipolar ionized nebula G316.80-00.05 oriented
approximately northeast to southwest as well an older outflow extending to the
east (Dalgleish et al., 2018). The parent molecular cloud in this region is
primarily in the form of a ridge that hosts the embedded massive stellar
cluster associated with the H ii region. The ridge is seen in absorption in
the visible, near-infrared, and at Spitzer-GLIMPSE wavelengths (Fig. B.1a in
Samal et al., 2018), and is clearly detected in emission at the Herschel-SPIRE
wavelengths and in the 870 µm ATLASGAL images. In this system the S110 bubble
is the northern lobe, which extends perpendicular to the filament, while S109
and S111 form the southern lobe (Fig. 1). The complex morphology of the region
is attributed to the inhomogeneous density distribution that causes the
ionization front to be distorted and enables leakage of ionizing photons to
large distances from the embedded cluster as evidenced by the presence of
diffuse emission from the ionized gas in S110 and S111 (Samal et al., 2018).
The dust temperature distribution derived from the Herschel Hi-GAL (Molinari
et al., 2010) continuum observations suggests that while the parental filament
is cold ($\sim 17$ K), going down to 14 K, the photon dominated regions (PDRs)
register temperatures well above 20 K (Samal et al., 2018). The column density
at the position of peak temperature position is $\sim 10^{23}$ cm-2, which
corresponds to a visual extinction of $\sim 100$ mag, which explains why the
stellar cluster ionizing the region is not detected in the visible and near-
IR.

In this paper, we refer to the ridge as G316.75 and the H ii region as
G316.80-00.05. The G316.75 ridge located at a distance of 2.69$\pm$0.45 kpc
(Watkins et al., 2019) extends over 13.6 parsec. The southern part of the
ridge, which is associated with an infrared dark cloud (SDC316.786-0.044
Peretto & Fuller, 2009), is more active than the northern part (Watkins et
al., 2019). Shaver et al. (1981) were the first to confirm the near kinematic
distance of G316.75 and determined that the source of the H ii region to be an
O6-type star. Using radiative transfer modeling of near- to far-infrared dust
continuum emission, Vig et al. (2007) had estimated that two O-type stars of
masses 45 $M_{\odot}$ and 25 $M_{\odot}$ are responsible for most of the
infrared luminosity of G316.75. This is also corroborated by high-angular-
resolution radio continuum and hydrogen recombination line observations
(Longmore et al., 2009) which detected two 24 GHz continuum sources at
resolutions better than 10″. Watkins et al. (2019) have subsequently concluded
that two high-mass stars with spectral types of O8.5–O7 V and O6.5–O6 V can
explain the observed emission from the ionized gas from the region. Dalgleish
et al. (2018) have studied the kinematics of the ionized gas using radio
recombination line (RRL) emission, and concluded that the strong velocity
gradient they observe could be a relic of the cloud’s initial angular
momentum. Study of the southern part of the ridge G316.75 is rendered
interesting by the fact that is a very active and young star-forming region,
harboring water, hydroxyl, and methanol masers, several O-stars, and a compact
X-ray source, and it is characterized by very dynamic gas conditions.

Figure 1: Bipolar nebula G316.80-00.05. (a) Composite color image of the
G316.790-0.045 complex: red, green, and blue indicate the Spitzer 24, 8.0, and
4.5 µm data, respectively (the 24 µm emission is saturated in the center). (b)
Composite color image of the complex: red is for the Herschel 250 µm emission
showing the cold filament at the waist of the nebula, green is for the 8 µm
PAH emission, and blue is for the 4.5 µm stellar emission. The three bubbles
listed in the (Churchwell et al., 2006) catalog are marked. Figure reproduced
from Samal et al. (2018) with permission.

Beyond the ionized gas, a significant fraction of the neutral (atomic and
molecular) material in the vicinity of such embedded clusters in the bubbles
is impacted by the far-ultraviolet (FUV:6–13 eV) radiation creating the so-
called PDRs. In the far-infrared (FIR) the most important cooling lines from
PDRs are the fine structure lines of [C ii] at 158 µm and [O i] at 63 and 145
µm, and to a lesser extent, high-$J$ CO lines, while polycyclic aromatic
hydrocarbon (PAH) emission and H2 lines dominate in the near- and mid-IR. The
ionization potential of carbon being 11.26 eV, C+ is expected to exist in
ionized, atomic, and molecular gas. As a result [C ii] has a unique position
as a tracer of PDR gas owing to the ubiquity of carbon and the moderate
excitation conditions (Eup=92 K, $n_{\rm cr}$=3000 cm-3) of the transition.
Thus, [C ii] at 158 µm is an ideal tool for tomographic study of the
morphology and kinematics of photoirradiated bubbles. The availability of far-
infrared observing facilities with multi-pixel receivers has recently enabled
study of the velocity fields in the Galactic PDRs using the information
derived from mapping observations of [C ii] at 158 µm.

In this paper we study the morphology and kinematics of the PDR gas in the
G316.75 region that is exposed to the FUV radiation from the embedded massive
stellar cluster that also gives rise to the bipolar H ii nebula G316.8-0.05.
For this purpose we primarily use the velocity information derived from
observations of [C ii] emission at 158 µm in combination with 13CO(2–1) and
HCO+(1–0) emission that arises primarily from the molecular material.

##  2 Observations

The paper makes use of several sets of publicly available data, which are
described below.

###  2.1 GREAT/SOFIA observations

We have retrieved observations of the 2P3/2$\rightarrow$2P1/2 fine structure
transition of ionized carbon (C+) at 1900.5369 GHz (157.74 µm) of the region
around the G316.8-00.05 from the data archive of the Stratospheric Observatory
for Infrared Astronomy (SOFIA; Young et al., 2012). The observations (Id:
06_0074; PI: J. Jackson) were carried out using the German REceiver for
Astronomy at Terahertz frequencies (GREAT; Heyminck et al., 2012) on 2018 June
22. The beam size for [C ii] was 14$\farcs$1\. The [C ii] map extends over a
500″$\times$600″ area and observations were carried out in total power on-the-
fly mode. We used the Level 4 data available on the SOFIA archive for the
region. All data were smoothed to a spectral resolution of 0.25 km s-1
resulting in an rms of 0.8 K. The [O i] 63 µm fine structure transition was
also mapped simultaneously as part of the same project. However, since the
mapping was optimized for [C ii], the [O i] data were undersampled and smeared
in the scan direction and hence were not used in this work.

###  2.2 SEDIGISM

We have used the data cubes for $J$=2–1 transition of 13CO and C18O that were
observed as a part of the large-scale (84 deg2) spectroscopic survey of the
inner Galactic disk, named Structure, Excitation and Dynamics of the Inner
Galactic Interstellar Medium (SEDIGISM; Schuller et al., 2021). These data
were observed with the APEX telescope between 2013 and 2016 with an angular
resolution of 30″ and a 1$\sigma$ sensitivity less than 1.0 K at 0.25 km s-1
velocity resolution.

###  2.3 MALT90

The Millimetre Astronomy Legacy Team 90 GHz (MALT90) survey (Foster et al.,
2011, 2013; Jackson et al., 2013) mapped a sample of high-density cores in 16
transitional lines including the high-density tracer HCO+(1–0) at 89188.5234
MHz. These data were taken with the 22 m Mopra radio telescope by individually
observing each of the $\sim 2000$ targeted clumps in a 3′$\times$3′ data cube.
The data cubes have an angular resolution of 38″ and an rms noise of 0.2 K at
a spectral resolution of $\sim 0.11$ km s-1. Three MALT90 observations lie
within the region studied here. We have in particular used the data for
HCO+(1–0) with a critical density of 2$\times 10^{5}$ cm-3 to understand the
kinematics of the underlying dense gas in the region.

##  3 Impact of Massive O-stars in G316.75 on the Ridge

Watkins et al. (2019) have performed a detailed analysis of the radio
continuum emission from the embedded G316.75 stellar cluster that irradiates
the ridge and its vicinity. The mass of the stellar cluster is 930$\pm$230
$M_{\odot}$ and it consists of four O-type stars one of which has a mass
larger than 48 $M_{\odot}$. The ionizing stars responsible for the two radio
continuum peaks seen in the Sydney University Molongolo Sky Survey (SUMSS
Mauch et al., 2003) were further constrained to be due to two high-mass stars
with spectral types O8.5–O7 and O6.5–O6. The dynamical age of a G316.75-like H
ii region with a radius of 6.5 pc and a Lyman continuum photon production
rate, $N_{\rm Lyc}$ = 1049.63 s-1 was estimated to be $\sim 2$ Myr.

Earlier studies using primarily tracers of ionized and molecular material
concluded that the G316.75 ridge is mostly unaffected by ionization owing to
the large electron and H2 gas densities in the immediate vicinity of the
ionizing high-mass stars (Watkins et al., 2019). These authors argued that
ionization due to OB stars happened after the formation of the molecular
ridge. Here we concentrate on the distribution and velocities of the FUV-
irradiated gas as traced by the PDR tracer – the [C ii] line at 158 µm – to
study the radiative and mechanical feedback of the ionizing stars primarily on
the S111 bubble.

##  4 Results

Figure 2: Comparison of [C ii] emission at 158 µm emission integrated between
-65 and -20 km s-1 (contours) with emission due to other tracers of dust,
PAHs, and molecular gas (in color and names marked in each panel). The [C ii]
contours are drawn from 90 to 540 K km s-1 in steps of 50 km s-1. The color
scale corresponding to each panel is included to the right of it. The units
are MJy sr-1 for the 8 µm image, Jy pixel-1 for 3″ pixels for the PACS 70 µm
image, and Jy beam-1 with a beam size of 21 ″ for the ATLASGAL 870 µm. The
color scales of spectroscopic images are all in units of K km s-1. The
velocity intervals of integration of 13CO(2–1), C18O(2–1) and HCO+(1–0) are
-60 to -10 km s-1, -47 to -20 km s-1 and -48 to -26 km s-1, respectively. For
each of the tracers shown as color images the beams are shown in the top left
corner of each panel. The radio continuum sources are marked with ”+” (white),
the ATLASGAL clumps (Csengeri et al., 2014) are marked with filled triangles
(black), the 70 µm bright far-infrared sources are marked with ”+” (black) and
the 70 µm dark FIR sources are marked with asterisks. The positional offsets
in arcseconds on the two axes are relative to the center ($\alpha_{2000}$:
14h45m20.4s, $\delta_{2000}$: -59°52′03$\farcs$5).

The velocity-integrated [C ii] map shows three bright emission regions, one
centered on the bright radio continuum emission in the southern part of the
ridge, one to the southeast of it, and one to the east, skirting the S111
bubble (Fig. 2). All three emission regions are connected by faint [C ii]
emission arising from diffuse PDR gas. When compared with the 8 µm Spitzer
image the [C ii] emission is clearly seen to delineate the rims skirting the
bubbles S109 and S111 in particular and matches well with the emission
features in the 8 µm which are due to FUV-excited PAH molecules in the PDRs.
The [C ii] emission also matches well with the far-infrared dust continuum
emission at 70 µm. The tracers of cold dust (870 µm ATLASGAL) and ambient
molecular material (13CO(2–1) and C18O(2–1)) primarily show the molecular
ridge. The southeastern [C ii] peak is detected in both 13CO and C18O,
suggesting a high column density of gas in the region. The HCO+(1–0) map
tracing densities in excess of 105 cm-3 shows that this line is detected in
the three prominently [C ii]-emitting regions, and the emission continues to
increase to the north toward the denser part of the molecular ridge. The
region hosts a large number of far-infrared sources, some of which are dark
even at 70 µm, and there are quite a few additional cold clumps detected in
the 870 µm ATLASGAL data (marked on Fig.2). The detection of HCO+(1–0)
particularly toward the easternmost parts of the northern rim surrounding S111
and the detection of associated far-infrared sources are clear indications of
ongoing star formation in dense clumps.

Figure 3: Maps of [C ii] emission (contours) in 1 km s-1 wide channels
centered at velocities indicated in the panels overplotted on the Spitzer 8 µm
continuum image. Contours (red) in the bottom two panels show the 843 MHz
radio continuum observed by SUMSS and are drawn from 0.2 to 4.2 Jy beam-1 in
steps of 0.5 Jy beam-1 (for a beam size of 52″$\times$45″).

##  5 Velocity distribution of [C ii]-emitting gas

###  5.1 Velocity-Channel Maps

Based on velocity-channel maps (Fig. A.1) we identify four distinct velocity
ranges of the [C ii] emission. Figure 3 compares the [C ii] emission in these
four velocity ranges centered at -45, -39, -33 and -24 km s-1 with the 8 µm
continuum image, which reveals the structure of the region very well. The
spatial segregation of emission in the four velocity regions indicates the
relative motion primarily of the rims of the bubble S111 and to some extent
that of S109. We find that [C ii] shows multiple velocity components spanning
an interval of 20 km s-1 spread over about 300″. As detailed later in the
paper, these velocity components can be explained by the expansion of the H ii
region S111.

###  5.2 Position–Velocity Diagrams

Figure 4: Velocity-integrated map of [C ii] emission with intensities ranging
between 10 and 540 km s-1 (shown by both color and contours). The contours are
separated by 50 K km s-1. Straight lines (red) with arrows (marked 1, 2, and
3) show the directions along which position-velocity diagrams shown in Fig. 5
are derived. Clouds A, B, and C discussed in the text are also marked. Also
shown is an ellipse approximating the S111 bubble traced by the [C ii]
contours. The ellipse has semi-major and semi-minor axes of 200″ and 125″
respectively and a position angle of 45° and is centered at the position with
an offset -100″,60″.

A more detailed view of the kinematics of the region can be obtained by
studying the position-velocity (p-v) diagrams along a few selected directions.
For ease of reference we identify the [C ii] prominent emission peaks as
clouds A, B, and C (Fig. 4). We considered directions along the ridge (Cut 1),
across the S111 bubble (Cut 2) and direction perpendicular to the ridge
connecting the [C ii] peaks (Cut 3) (Fig.4). We compare the p-v diagrams of
13CO(2–1) and [C ii] along the three directions. The 13CO(2–1) emission is
representative of the ambient molecular material, while the [C ii] mostly
represents FUV-irradiated neutral gas except for the immediate neighborhood of
the embedded cluster G316.8-00.05 (Fig. 5). The dotted line in all the p-v
diagrams denotes the centroid velocity of the molecular material, supported
also by the C18O(2–1) p-v plots (Fig. A.2). Along Cut 1, in the region near
the embedded cluster the [C ii] emission compares fairly well at velocities
around -39 km s-1 but shows a distinctly red-shifted component reminiscent of
an expanding bubble. The same feature also shows up as isolated red-shifted
blobs on the northeastern end of Cut 3. We identify this as the interface of
the bubbles S109/S111 that is expanding away from the observer and propose
that the expansion is restricted in the opposite direction by the high density
molecular material that forms the ridge. Considering the difference in
velocities of the molecular material and the PDR gas traced by the [C ii], we
estimate an expansion velocity of around 7 km s-1 for part of the rim detected
around S111. Along Cut 2 across the bubble S111, the northern [C ii] emission
feature is faint and compact in velocity in 13CO(2–1). However, [C ii] at the
same position is significantly broadened around the centroid velocity of the
ambient medium, possibly due to the presence of young stellar objects (YSOs)
that are very bright in the infrared. Along Cut 2, at the inner edge of the
bubble S111, the [C ii] emission is markedly red-shifted. Along Cut 3, the
southwestern clump shows red-shifted emission, while the northeastern clump is
somewhat blue-shifted relative to the velocity of the ambient molecular cloud.
The velocity distribution seen in the p-v diagrams of [C ii] is consistent
with emission from FUV-irradiated gas at velocities suggestive of an expanding
bubble.

Figure 5: Position–velocity maps of 13CO(1–0) and [C ii] emission along the
Cuts 1, 2, and 3 marked in Fig.4.

###  5.3 Velocity distribution seen in Azimuthally averaged spectra

Figure 6: Plot of velocity distribution as a function of radial distance from
the position (-100″,60″). Each line of the map is the radial profile for one
particular plane of the cube and is created by computing annular averages on
the data cube. The horizontal dashed lines are drawn at distances of 50″,
100″, 135″, 180″, and 250″ from (-100″,-60″) and show the radial distances of
the clouds A, B and C from the same point.

The analysis of [C ii] emission as a function of velocity in the previous
sections of the paper suggests that [C ii] traces parts of the expanding
bubble S111. The shape of the bubble seen in [C ii] is an ellipse with a size
of 400″$\times$250″centered at $\sim$ -100″, +60″, with a position angle of
45°(Fig 4). In order to improve the signal-to-noise ratio of the high-velocity
[C ii] emission, we have computed the [C ii] emission as a function of the
radial distance from the assumed center of the bubble, i.e. at an offset
(-100″,60″) (Fig. 6). Each line of the map is the radial profile for one
particular plane of the cube and is created by computing annular averages on
the data cube. This plot shows that the [C ii] emission is quite broad over
the whole nebula, with line widths in the range of 18 – 22 km s-1 (Fig. 6)
consistent with an expanding lobe/bubble with a radius of approximately 220″.
The [C ii] emission is particularly strong in the interface regions of clouds
A, B, and C. There is a clear velocity gradient between clouds A and C, with
the bulk of the gas being more redshifted in cloud A than in clouds B and C,
likely due to the interaction between the S111 and S109 bubbles. In these
radially averaged data the broadest [C ii] lines are seen in the direction of
cloud C, whereas the p-v plots (Fig. 5) show that the broadest [C ii] emission
is in the direction of cloud B, which hosts the O-stars responsible for
ionizing the H ii region. This can be explained by the fact that the cloud C
is much more extended in [C ii] than cloud B, so that it dominates the radial
average of the data, which is more sensitive to the extended emission.

###  5.4 Analysis of [C ii] spectra at selected positions

Figure 7 shows a comparison of [C ii], 13CO(2–1) and C18O(2–1) spectra at
selected positions where [C ii] emission peaks are observed and which coincide
with the locations of Hi-GAL far-infrared sources. The selected positions also
sample the environment in the different parts of the arc-shaped rim around
S111. The [C ii] spectra trace the warm gas around the embedded FIR sources as
well as the gas photodissociated by the embedded cluster G316.8-0.05.
Therefore the emission profiles appear more complicated than the 13CO and C18O
profiles, which are dominated by the dense molecular gas. Table 1 presents the
results of single- or multiple-component Gaussian fits to the [C ii],
13CO(2–1) and C18O(2–1). We find that the C18O(2–1) profiles are well fitted
by a single velocity component although the central velocity appears to vary
by more than 5 km s-1 across the bubble. At all positions, one of the
components seen in [C ii] matches with the component seen in C18O and with one
of the two velocity components typically seen in 13CO. The larger line widths
(5–12 km s-1) seen in [C ii] are likely due to photoevaporative flows of
possibly overlapping PDRs. The broad lines seen in 13CO(2–1) and HCO+(1–0)
likely arise from outflows associated with the young embedded far-infrared
sources. In particular, the outflow from the embedded YSO is clearly detected
in the 13CO and HCO+ spectra of the position with an offset of (-20″, 10″)
which is located in cloud A (Fig. 7).

The [C ii] spectra at all positions except toward the center of the star
cluster at (55″, 155″) are more redshifted than the molecular emission. The
hydrogen RRLs observed by Longmore et al. (2009) have central velocities
between -42 and -46 km s-1, i.e., they are significantly blue-shifted compared
to the [C ii], 13CO, and C18O spectra. However, the RRL spectra are much
broader than the molecular lines and [C ii]. The line widths are between 20
and 24 km s-1, and therefore they overlap with some of the [C ii] spectra.
While [C ii] may originate from ionized gas as well, in all cases where [C ii]
is detected at blue-shifted velocities, 13CO(2–1) is also detected, which
implies that the [C ii] emission is primarily associated with the neutral gas.

Figure 7: Comparison of spectra at selected positions in the observed [C ii]
map (shown in the middle). In each panel the offset of the position relative
to ($\alpha_{2000}$: 14h45m20.4s, $\delta_{2000}$: -59°52′03$\farcs$5) is
marked and [C ii] (black), CO(1–0) (purple), 13CO(2–1) (red), C18O(2–1)
(blue), and HCO+(1–0) (brown) spectra are shown. Single- or multiple-component
Gaussian fits to the three spectra at each position were performed and the
results are shown as green continuous curves with details presented in Table
1. The spectra are scaled by the factors mentioned in the panels for better
visibility.  Table 1: Results of Gaussian fits to spectra

Position | Transition | $\upsilon_{\rm LSR}$ | $\Delta\upsilon_{\rm LSR}$ | $I$ |  $T_{\rm kin}$a
---|---|---|---|---|---
(arcsec) |  | (km s-1) | (km s-1) | (Kkm s-1) | (K)
(20, 100) | [C ii] | -41.52$\pm$0.07 | 5.07$\pm$0.21 | 59.97$\pm 1.89$ | 42
|  | -36.04$\pm$0.14 | 2.98$\pm$0.36 | 27.21$\pm 1.47$ |
|  | -23.88$\pm$0.15 | 5.55$\pm$0.35 | 12.70$\pm 1.50$ |
|  13CO(2–1) | -38.92$\pm$0.49 | 6.90$\pm$0.28 | 25.37$\pm 1.42$ |
|  | -40.03$\pm$0.18 | 6.77$\pm$0.20 | 66.14$\pm 1.55$ |
| C18O(2–1) | -39.40$\pm$0.24 | 6.71$\pm$0.43 | 18.25$\pm 1.22$ |
(55, 155) | [C ii] | -39.91$\pm$0.25 | 6.85$\pm$0.25 | 281.08$\pm 3.84$ | 90
|  | -40.11$\pm$0.25 | 3.30$\pm$0.25 | 47.03$\pm 3.84$ |
|  | -22.55$\pm$0.25 | 5.43$\pm$0.25 | 35.76$\pm 3.84$ |
|  13CO(2–1) | -38.69$\pm$0.5 | 3.64$\pm$0.5 | 39.88$\pm 2.93$ |
|  | -39.57$\pm$0.5 | 6.83$\pm$0.5 | 73.22$\pm 2.93$ |
| C18O(2–1) | -38.05$\pm$0.08 | 4.40$\pm$0.19 | 31.08$\pm 1.08$ |
(-20, 10) | [C ii] | -34.91$\pm$0.03 | 2.38$\pm$0.08 | 47.51$\pm 2.39$ | 88
|  | -33.94$\pm$0.03 | 8.39$\pm$0.09 | 271.40$\pm 2.75$ |
|  13CO(2–1) | -35.41$\pm$0.04 | 3.16$\pm$0.12 | 38.82$\pm 1.80$ |
|  | -38.04$\pm$0.31 | 13.75$\pm$0.60 | 48.01$\pm 2.47$ |
| C18O(2–1) | -35.55$\pm$0.14 | 2.96$\pm$0.35 | 8.98$\pm 0.88$ |
(-140, 185) | [C ii] | -38.98$\pm$0.02 | 3.63$\pm$0.08 | 144.18$\pm 6.36$ | 93
|  | -34.71$\pm$0.09 | 8.02$\pm$0.13 | 297.54$\pm 6.69$ |
|  | -42.64$\pm$0.02 | 2.67$\pm$0.05 | 67.73$\pm 1.43$ |
|  13CO(2–1) | -36.31$\pm$0.21 | 4.55$\pm$0.48 | 17.98$\pm 1.83$ |
|  | -40.70$\pm$0.03 | 3.22$\pm$0.07 | 69.32$\pm 1.71$ |
| C18O(2–1) | -40.71$\pm$0.06 | 2.11$\pm$0.15 | 10.15$\pm 0.59$ |
(-130, 85) | [C ii] | -34.09$\pm$0.25 | 4.37$\pm$0.25 | 8.03$\pm 2.56$ | 69
|  | -34.03$\pm$0.25 | 9.22$\pm$0.25 | 195.84$\pm 2.56$ |
|  | -32.59$\pm$0.25 | 3.34$\pm$0.25 | 41.51$\pm 2.56$ |
|  13CO(2–1) | -32.51$\pm$0.06 | 3.28$\pm$0.13 | 26.00$\pm 0.92$ |
|  | -39.55$\pm$0.13 | 4.72$\pm$0.39 | 21.04$\pm 1.19$ |
| C18O(2–1) | -32.58$\pm$0.18 | 2.50$\pm$0.37 | 5.25$\pm 0.71$ |

Note. — a Planck-corrected peak temperature of the [C ii] spectrum assumed to
be optically thick.

##  6 Column Density and Mass of the [C ii]-emitting gas

The velocity distribution of the region with the strongly red-shifted
component detected in the [C ii] emission does not allow for the detection of
the much fainter hyperfine component of [13C ii] line, which lies at
$\upsilon=+11$ km s-1 relative to the [C ii] line. However, in recent
velocity-resolved observations, the [C ii] emission from most Galactic PDRs
has been found to be optically thick (Mookerjea et al., 2021; Guevara et al.,
2020; Mookerjea et al., 2019). The PDR gas in the G316.75 region being studied
here is reasonably bright and high-density clumps are seen in the [C ii]
emission as well. We thus consider the [C ii] emission to be optically thick
and use the Rayleigh-Jeans corrected peak brightness temperature of [C ii]
spectra as an estimate of the kinetic temperature of the PDR gas (Column 6 of
Table 1) and also for the estimation of $N$(C+) from the observed velocity-
integrated [C ii] intensity. We first consider the contribution of ionized gas
close to the H ii region. Comparison of the [C ii] emission with the
H57$\alpha$ (Longmore et al., 2009) emission map suggests that the peak of the
[C ii] emission is shifted by about 10″, which is within the 14″ beam of the
[C ii] data. Watkins et al. (2019) explained the small ionized gas fraction in
the ridge G316.75 as as a consequence of the large electron ($n_{e}$=5000
cm-3) and H2 gas densities in the immediate vicinity of the ionizing high-mass
stars. If we assume $n_{e}$=5000 cm-3 and follow the formulation presented by
Pabst et al. (2017, Equation 7) to estimate the contribution of the ionized
gas to the [C ii] intensities close to the H ii region we obtain an integrated
[C ii] intensity of 497 K km s-1. The [C ii] intensity estimated to arise from
the ionized gas is therefore comparable to the peak [C ii] intensity of 437 K
km s-1 found at this position. However, unlike the prediction of the model
constructed by Watkins et al. (2019), the H57$\alpha$ emission suggests that
such high-density ionized gas to be confined to within 30″ of the peak of
ionized emission. At the other extreme ($n_{e}\sim 60$ cm-3) of the range of
electron densities proposed by Watkins et al., the [C ii] emission would be
only 6 K km s-1, far less than the faintest level (180 K km s-1) level of [C
ii] emission in the immediate vicinity of the H ii region (Fig. 4) . Thus, we
conclude that while for the peak position of the H ii region there is a
significant contribution of the ionized gas to the [C ii] emission, the more
extended emission is dominated by the neutral PDR gas.

In the neutral region, where C+–H and C+–H2 collisions dominate, we estimate
$N$(C+) (in cm-2) from the observed integrated [C ii] intensities using Eq.
(26) from Goldsmith et al. (2012) as follows:

$N(C^{+})={\rm 2.91\times 10^{15}\left[1+0.5{\rm
e^{91.25/T_{kin}}}\left(1+\frac{A_{ul}}{C_{ul}}\right)\right]\int T_{\rm
mb}d\upsilon}$ (1)

where $A_{\rm ul}$ = 2.3$\times 10^{-6}$ s-1, $T_{\rm kin}$ is the gas kinetic
temperature, the collision rate is $C_{\rm ul}$ = $R_{\rm ul}n$ with $R_{\rm
ul}$ being the collision rate coefficient with H2 or H0, which depends on
$T_{\rm kin}$, and $n$ is the volume density of H. At $T_{\rm kin}$ = 100 K
for C+–H collisions $R_{\rm ul}$ = 7.6$\times 10^{-10}$ cm3 s-1 and $R_{ul}$ =
3.8$\times 10^{-10}$ cm3 s-1 for C+–H2 collisions. Since the critical density
of the [C ii] transition is $n_{\rm cr}$=3000 cm-3, and since it is likely
that most of the [C ii] detected could be at such densities along with some
emission arising from clumps with densities exceeding 105 cm-3, we assume a
density of 104 cm-3 to estimate the $N$(C+) density distribution of the
region. Assuming a kinetic temperature of 100 K for C+ and a density of
$n$=104 cm-3, and excitation from C+–H2 collisions, the column density is
found to lie in the range (0.7–4.5)$\times 10^{18}$ cm-2. A combination of
slightly different $T_{\rm kin}$ and density and the inclusion of C+–H
collisions would result in an uncertainty of no more than 30% in the column
density. Based on the values of $A_{\rm ul}$ = 2.3$\times 10^{-6}$ s-1 and
$n_{\rm cr}$ at 100 K and C/H = 3$\times 10^{-4}$, the [C ii] line is expected
to be optically thick for $N$(C+) = 4$\times 10^{17}$ cm-2 (Table 2.7 in
Tielens, 2005). The estimated values of $N$(C+) thus validate our assumption
of the [C ii] line being optically thick.

We estimate the mass of the [C ii]-emitting gas, which is dominated by the
three regions named clouds A, B and C (Fig.4). For the observed [C ii]
intensities, a distance of 2.7 kpc and a C/H abundance ratio of 1.3$\times
10^{-4}$ of which approximately 50% of carbon resides in C+, we find the total
molecular H2 masses of clouds A, B and C to be 1000, 640, and 2500 $M_{\odot}$
respectively. We have also estimated the mass of the molecular gas in clouds
A, B, and C using the $N$(H2) maps derived by Samal et al. (2018) based on
Herschel observations. We find the H2 mass estimated from dust continuum
observations in clouds A, B, and C to be 2000, 3400, and 5000 $M_{\odot}$
respectively. The estimate of the total molecular gas mass within the lowest
[C ii] contour shown in Fig 4 is 6900 $M_{\odot}$ from [C ii] data and 18,100
$M_{\odot}$ from the dust continuum. Thus, we see that in clouds A and C,
which are bright in far-infrared dust continuum as well as [C ii], about 50%
of the mass is traced by the [C ii] emission, whereas in cloud B, in which the
ionizing stars reside, only 18% of the mass is traced by [C ii]. The lower
fraction of molecular mass being traced by [C ii] in this region is likely due
to the high densities in the ridge, which shields it from the FUV radiation,
unlike the rest of the cloud. The [C ii] emission seen from the ridge is
likely due to dense the ionized gas near the center of the H ii region and
diffuse PDR further out. This indication of higher density material to the
north along the ridge is corroborated by the analysis presented by Watkins et
al. (2019) in support of the lower than expected effect of the embedded
cluster on the ridge itself. We propose that the bipolar H ii region clears
off part of the bubble, and the FUV radiation leaking through the bubble
creates PDRs on the surface of the G316.75 ridge.

##  7 Estimate of Mechanical feedback from the stars on S111

###  7.1 Energy budget in S111

The triple bubble system of S109, S110 and S111 has putatively been created by
the massive stellar cluster embedded in the G316.75 ridge. While the three
bubbles are connected through their history of formation, in the present
configuration the H ii regions appear to be segregated by the dense ridges or
shells of swept-up gas and hence can be studied in isolation. Our study
focuses on the mechanical and radiative feedback from the embedded stars.

We first estimate the role of mechanical feedback from the ionizing stars by
comparing the kinetic energy of the shell to the the thermal energy of the
ionized gas in S111. The [C ii] observations suggest that the denser parts of
the clouds A and C represent part of the shell created by the expanding H ii
region. These have velocities red-shifted by about 5–6 km s-1 relative to
cloud B, which hosts the ionizing stars (Table 1, Fig.5). Additionally more
diffuse PDR material also traced by [C ii] are red-shifted by 15 km s-1 at the
boundary between the bubbles S109 and S111. Thus, the shells of S111 appear to
move at relatively slow speeds, particularly when contrasted with the shells
that have been observed in Orion, RCW 120, and RCW 49 (Pabst et al., 2020;
Luisi et al., 2021; Tiwari et al., 2021). The high-density material toward the
northern part of the ridge seen in the molecular line observations most likely
contributes to the slowing down of the expansion of the H ii region and the
shell. We estimate the kinetic energy of the clouds A and C with a mass of
7000 $M_{\odot}$ moving with a velocity of 7 km s-1 relative to B to be
3.4$\times 10^{48}$ erg.

The H ii region S111 has a radius of $\sim$ 160″ (2 pc at a distance of 2.7
kpc) (Fig. 4) and Watkins et al. (2019) estimated a temperature $T_{\rm HII}$
= 6600 K and density $n_{\rm e}$=60–5000 cm-3. The total thermal energy of the
S111 H ii region can be written as ${\rm
E_{th}=\frac{3}{2}n_{e}\frac{4\pi}{3}\,R^{3}k_{B}T_{\rm HII}}$, where ${R}$ is
the radius of the H ii region and $k_{\rm B}$ is the Boltzmann constant. We
estimate the $E_{\rm th}$ to range between 8.1$\times 10^{46}$ and 6.7$\times
10^{48}$ erg for $n_{\rm e}$=60 and 5000 cm-3 respectively.

Thus, depending on the actual electron density of the H ii region, the kinetic
energy of the molecular shells primarily detected in [C ii] is between $\sim$
0.5–40 times the thermal energy of the H ii region.

###  7.2 Pressure balance of S111

The impact of high-mass stars on the ambient molecular cloud can be further
quantified in terms of the pressures associated with direct radiation ($P_{\rm
dir}$), dust-processed radiation ($P_{\rm IR}$), photoionization heating
($P_{\rm ion}$) and shock heating from stellar winds ($P_{\rm X}$). No X-ray
mapping observations are available for the H ii region associated with S111,
hence $P_{\rm X}$ could not be estimated. Thus, we next consider the pressures
associated with the other feedback mechanisms from the high-mass stars in the
neighborhood of S111.

The pressure due to the photoionization heating of the gas in the region,
$P_{\rm ionized}/k_{\rm B}=n_{\rm e}T_{\rm HII}$ is estimated to be between
(0.4–33)$\times 10^{6}$ K cm-3 for electron densities of 60–5000 cm-3
respectively.

Following Olivier et al. (2021), we define the direct radiation pressure
($P_{\rm dir}$) as the momentum available to drive motion in the shells in the
H ii region at a radius $R$ from the central stars (Eq. 3)

${\rm P_{dir}}={\rm\frac{3L_{bol}}{4\pi R^{2}c}},$ (2)

where $L_{\rm bol}$ is the bolometric luminosity of central stars. The total
bolometric luminosity of the two central ionizing stars, one O7V and the other
O6V, is 3.5$\times 10^{5}$ $\mathrm{L}_{\odot}$ and the total production rate
of ionizing photons is 1.56$\times 10^{49}$ s-1. The direct pressure ($P_{\rm
dir}/k_{\rm B}$) at the location of the inner edge of the shell located at a
radius of 160″ (2 pc) from the stellar cluster is therefore 9.3$\times 10^{5}$
K cm-3.

In another feedback mode, stellar radiation is absorbed by dust and thermally
reradiated in the IR wavelengths, and this induces a dust-processed radiation
pressure that is given by (Eq. 2 in Olivier et al., 2021)

${\rm P_{IR}}=\frac{1}{\rm 3V}a{\rm T_{d}^{4}}\frac{4\pi}{3}{\rm(R_{\rm
shell}^{3}-R^{3})},$ (3)

where $R_{\rm shell}$ is the outer radius of the shell, $R$ is the inner
radius of the shell, which is the same as the radius of the H ii region, and
$V$ is the volume of the H ii region. Based on the dust temperatures $T_{\rm
d}$ derived by Samal et al. (2018), we assume the dust temperature in the
shell to be around 24 K. From the position-velocity plots (Fig 5) we estimate
the thickness of the shell to be 20″ (0.3 pc), resulting in a dust-processed
radiation pressure $P_{\rm IR}/k_{\rm B}$ of 6$\times 10^{6}$ K km s-1.

Finally, we estimate the thermal pressure ($P_{\rm PDR}$) of [C ii]-emitting
PDR gas. We assume the hydrogen volume density of the PDR to be $n_{\rm H}$=
3$\times 10^{3}$ cm-3, which is the critical density ($n_{\rm cr}$) for the [C
ii] transition, and the gas temperature ($T_{\rm PDR}$) to be 100 K based on
our estimates in Section 6. Thus ${\rm P_{\rm PDR}}=n_{\rm H}k_{\rm B}T_{\rm
PDR}=3\times 10^{5}$ K km s-1. It is likely that the value of $P_{\rm PDR}$ is
an underestimate, since we have considered only the diffuse component of the
PDR and it could easily be a factor of three larger, in which case it would be
comparable to $P_{\rm dir}$.

We find that among all factors contributing to pressure, $P_{\rm ion}$ is the
most important, and depending on the electron density it is of similar in
magnitude to $P_{\rm IR}$ of the shell around the H ii region. In this
analysis, we have not considered the dust-processed radiation pressure inside
the H ii region since there is very faint or almost no dust continuum detected
from within the H ii region even in the SPIRE 250 µm maps. We find that
$P_{\rm dir}$ and $P_{\rm PDR}$ are lower than $P_{\rm IR}$ by factors of 6
and 20, respectively. Olivier et al. (2021) considered resolved compact H ii
regions and found that most of the sources are dominated by $P_{\rm IR}$, and
the median $P_{\rm dir}$ and $P_{\rm ion}$ are smaller than the median $P_{\rm
IR}$ by factors of 6 and 9 respectively. Based on the radial dependence of
pressure terms, these authors concluded that the H ii regions transition from
$P_{\rm IR}$-dominated to $P_{\rm ion}$-dominated at radii of $\sim 3$ pc. Our
conclusions for the S111 H ii region, which has a radius of about 2 pc with a
shell of about 0.3 pc thickness, are consistent with the results obtained by
Olivier et al. (2021).

##  8 Discussion

The velocity-resolved observations of [C ii] have enabled a tomographic study
of the expanding shell associated with the bubble S111 formed by the embedded
cluster G316.8-00.05. The expansion of the bubble at $\sim 7$ km s-1 is
primarily observed in the red-shifted velocities with no blue-shifted
counterparts. This suggests that the bubble is able to expand more freely in a
direction away from the observer, possibly due to the presence of higher
density material between the rim of S111 and the observer. The H ii region
associated with S111 plays a key role in providing the pressure for the
expansion of the shell in the ambient molecular medium. We estimate ionization
of hydrogen ($P_{\rm ion}$) and dust-reprocessed radiation ($P_{\rm IR}$) to
be the primary contributors to the pressures in the region. We find that the
kinetic energy of part of the molecular shell associated with clouds A and C
is 0.5–40 times the total thermal energy of the H ii region. We do not see any
evidence of injection of energy by stellar winds from the massive stars; the
moderate expansion velocity of the shell can be well explained by the smooth
expansion of the H ii region. This is in contrast to several recent studies of
the feedback of massive stars on their environment, which imply a significant
role of the stellar winds in the Orion veil shell, RCW 120, and RCW 49 (Pabst
et al., 2020; Luisi et al., 2021; Tiwari et al., 2021). There is a distinct
difference between the ionizing stars in the G316.75 region and regions
studied by the above-mentioned works. The ionizing stars in Orion, RCW 120,
and RCW 49 have already cleared the surrounding region, whereas the
G316.80-0.05 cluster is still deeply embedded in dense molecular material with
a line-of-sight visual extinction exceeding 100 mag. This suggests that the
impact of the stellar winds, if any, from G316.8-0.05 could have also been
contained within a very small region due to the high densities of the region.

The dense rim created by the expansion of the bubble has been studied in terms
of the three clouds A, B, and C, with cloud B being coincident with the ridge
G316.75, which is dark in the mid-infrared and hosts the massive stellar
cluster that gave rise to the triple bubble system. Clouds A and C are created
by the compression of molecular material by the expanding bubble and contain
embedded far-infrared sources that are indicative of formation of a new
generation of stars. The two [C ii] emission peaks in cloud B are associated
with far-infrared sources detected at $\lambda>70$ µm as well as with 870 µm
ATLASGAL sources. The [C ii] peak in the cloud C coincides with the ATLASGAL
source G316.768-00.026, which remains undetected at $\lambda<160$ µm,
consistent with a YSO at a very early stage of evolution. Most of the YSOs
identified based on mid-infrared colors (Samal et al., 2018) lie along the
northern part of the region, with one lying very close to the head of cloud A,
one close to cloud B, and two in cloud C. One of the reasons for the
nondetection of YSOs in the region is attributed to the high degree of
obscuration along this line of sight.

Cloud A in particular appears to have a structure similar to the pillars
created by the UV erosion of molecular material (Fig. 3). At the head of cloud
A (-20″,10″) an outflow is detected in the form of a pronounced blue wing in
the 13CO(2–1) and HCO+(1–0) spectra (Fig. 7) arising from the embedded YSO.
The outflow coincides with the ATLASGAL source G316.778-00.096, which is also
detected in the far-infrared Hi-GAL images with its intensities peaking at 160
µm. Elia et al. (2017) have identified this source as HIGALBM 316.7799-0.0942
and have derived a mass of 234.6 $M_{\odot}$, a bolometric luminosity of 1187
$\mathrm{L}_{\odot}$, and a bolometric temperature of 38.6 K. While no other
signs of star formation have yet been reported, the detection of broad spectra
consistent with an outflow, high densities (based on detection of HCO+(1–0))
along with the presence of a compact far-infrared source are clearly
indicative of ongoing star formation activity. Based on the mass estimated for
the continuum source and the distance to the region, it is possible that a
star cluster is being formed at the head of cloud A. The H ii region of S111
is estimated to have an age of 2 Myr, which when compared with typical
formation timescales of high-mass stars implies that the star formation
activity in cloud A is more recent and hence likely triggered by the expansion
of S111.

##  9 Conclusions

We have used velocity-resolved observations of the
2P${}_{3/2}\rightarrow^{2}$P1/2 transition of [C ii] and the $J$=2–1
transition of 13CO and C18O to study the infrared bubble S111, which lies in
the region surrounding the southern part of the infrared dark ridge G316.75.
The [C ii] emission from the region primarily arises from the shells
irradiated by the embedded cluster G316.80–0.05, part of this photodissociated
gas in the shell in cloud C not being detected in 13CO(2–1) and C18O(2–1). In
the dense ridge G316.75 the [C ii] emission arises due to photodissociation of
the molecular material by the stellar radiation leaking out into the bubble
and not through the dense ridge. The [C ii] emission traces the low-density
PDR gas and shows dense clumps that emit HCO+(1–0) as well. However, in the
absence of any other PDR tracer it is not possible to further constrain the
properties of the PDR gas.

The velocity information of the [C ii] data in particular enabled us to
decipher the morphology of the expanding shell and evaluate the feedback from
the massive star cluster on the surrounding ISM. The expanding shell is
primarily propelled by the H ii region created by the embedded cluster
G316.8-0.05. Azimuthally averaged [C ii] data were used to conclusively show
that the entire velocity structure seen in [C ii] can be explained by layers
of a single shell expanding with a moderate velocity of $\sim 7$ km s-1 around
S111. We conclude that the kinetic energy of the expanding part of the shell
is comparable to or larger than the thermal energy of the H ii region. We have
also quantified the contribution of the stellar radiation to the pressure in
the region and conclude that the thermal pressure of the H ii region and the
dust-processed radiation pressure are the most significant components. We
present here a clear evidence of star formation triggered by the compression
of the ISM due to the expanding bubble S111.

The author acknowledges the insightful discussions with Dr G. Sandell during
the writing of the manuscript. This work was supported by funding by the
Department of Atomic Energy, Government of India, under Project Identification
No. RTI 4002. Based on observations made with the NASA/DLR Stratospheric
Observatory for Infrared Astronomy (SOFIA). SOFIA is jointly operated by the
Universities Space Research Association, Inc. (USRA), under NASA contract
NAS2-97001, and the Deutsches SOFIA Institut (DSI) under DLR contract 50 OK
0901 to the University of Stuttgart. The development of GREAT was financed by
the participating institutes, by the Federal Ministry of Economics and
Technology via the German Space Agency (DLR) under Grants 50 OK 1102, 50 OK
1103 and 50 OK 1104 and within the Collaborative Research Centre 956, sub-
projects D2 and D3, funded by the Deutsche Forschungsgemeinschaft (DFG). This
research has made use of the VizieR catalog access tool, CDS, Strasbourg,
France. The original description of the VizieR service was published in A&AS
143, 23. The ATLASGAL project is a collaboration between the Max-Planck-
Gesellschaft, the European Southern Observatory (ESO) and the Universidad de
Chile. It includes projects E-181.C-0885, E-078.F-9040(A), M-079.C-9501(A),
M-081.C-9501(A) plus Chilean data. This paper made use of information from the
SEDIGISM survey database located at https://sedigism.mpifr-
bonn.mpg.de/index.html, which was constructed by James Urquhart and hosted by
the Max Planck Institute for Radio Astronomy.

## References

  * Churchwell et al. (2006) Churchwell, E., Povich, M. S., Allen, D., et al. 2006, ApJ, 649, 759
  * Csengeri et al. (2014) Csengeri, T., Urquhart, J. S., Schuller, F., et al. 2014, A&A, 565, A75
  * Dalgleish et al. (2018) Dalgleish, H. S., Longmore, S. N., Peters, T., et al. 2018, MNRAS, 478, 3530
  * Deharveng et al. (2010) Deharveng, L., Schuller, F., Anderson, L. D., et al. 2010, A&A, 523, A6
  * Elia et al. (2017) Elia, D., Molinari, S., Schisano, E., et al. 2017, MNRAS, 471, 100
  * Foster et al. (2011) Foster, J. B., Jackson, J. M., Barnes, P. J., et al. 2011, ApJS, 197, 25
  * Foster et al. (2013) Foster, J. B., Rathborne, J. M., Sanhueza, P., et al. 2013, PASA, 30, e038
  * Goldsmith et al. (2012) Goldsmith, P. F., Langer, W. D., Pineda, J. L., & Velusamy, T. 2012, ApJS, 203, 13
  * Guevara et al. (2020) Guevara, C., Stutzki, J., Ossenkopf-Okada, V., et al. 2020, A&A, 636, A16
  * Hanaoka et al. (2020) Hanaoka, M., Kaneda, H., Suzuki, T., et al. 2020, PASJ, 72, 5
  * Heyminck et al. (2012) Heyminck, S., Graf, U. U., Güsten, R., et al. 2012, A&A, 542, L1
  * Jackson et al. (2013) Jackson, J. M., Rathborne, J. M., Foster, J. B., et al. 2013, PASA, 30, e057
  * Krumholz et al. (2018) Krumholz, M. R., Burkhart, B., Forbes, J. C., & Crocker, R. M. 2018, MNRAS, 477, 2716
  * Longmore et al. (2009) Longmore, S. N., Burton, M. G., Keto, E., Kurtz, S., & Walsh, A. J. 2009, MNRAS, 399, 861
  * Luisi et al. (2021) Luisi, M., Anderson, L. D., Schneider, N., et al. 2021, Science Advances, 7, eabe9511
  * Mauch et al. (2003) Mauch, T., Murphy, T., Buttery, H. J., et al. 2003, MNRAS, 342, 1117
  * Molinari et al. (2010) Molinari, S., Swinyard, B., Bally, J., et al. 2010, A&A, 518, L100
  * Mookerjea et al. (2019) Mookerjea, B., Sandell, G., Güsten, R., et al. 2019, A&A, 626, A131
  * Mookerjea et al. (2021) Mookerjea, B., Sandell, G., Veena, V. S., et al. 2021, A&A, 648, A40
  * Olivier et al. (2021) Olivier, G. M., Lopez, L. A., Rosen, A. L., et al. 2021, ApJ, 908, 68
  * Pabst et al. (2017) Pabst, C. H. M., Goicoechea, J. R., Teyssier, D., et al. 2017, A&A, 606, A29
  * Pabst et al. (2020) —. 2020, A&A, 639, A2
  * Peretto & Fuller (2009) Peretto, N., & Fuller, G. A. 2009, A&A, 505, 405
  * Samal et al. (2018) Samal, M. R., Deharveng, L., Zavagno, A., et al. 2018, A&A, 617, A67
  * Schuller et al. (2021) Schuller, F., Urquhart, J. S., Csengeri, T., et al. 2021, MNRAS, 500, 3064
  * Shaver et al. (1981) Shaver, P. A., Retallack, D. S., Wamsteker, W., & Danks, A. C. 1981, A&A, 102, 225
  * Tielens (2005) Tielens, A. G. G. M. 2005, The Physics and Chemistry of the Interstellar Medium (Cambridge University Press)
  * Tiwari et al. (2021) Tiwari, M., Karim, R., Pound, M. W., et al. 2021, ApJ, 914, 117
  * Vig et al. (2007) Vig, S., Ghosh, S. K., Ojha, D. K., & Verma, R. P. 2007, A&A, 463, 175
  * Watkins et al. (2019) Watkins, E. J., Peretto, N., Marsh, K., & Fuller, G. A. 2019, A&A, 628, A21
  * Young et al. (2012) Young, E. T., Becklin, E. E., Marcum, P. M., et al. 2012, ApJ, 749, L17

##  Appendix A Additional Figures

Figure A.1 presents the channel maps of [C ii] emission discussed in Section
5.1. Figure A.2 shows the position-velocity diagrams for C18O(2–1) which are
discussed in Section 5.2.

Figure A.1: Channel maps of [C ii] emission (contours) overlaid on Spitzer 8
µm continuum image (color). The positional offsets are relative to the center
($\alpha_{2000}$: 14h45m20.4s, $\delta_{2000}$: -59°52′03$\farcs$5).  Figure
A.2: Position-velocity maps of C18O(1–0) and [C ii] emission along the Cuts 1,
2 and 3 marked in Fig.4.