# The VLT MUSE NFM view of outflows and externally photoevaporating discs near
the Orion Bar††thanks: Based on observations collected at the European
Southern Observatory under ESO programme 110.259E.001.

Thomas J. Haworth1, Megan Reiter2, C. Robert O’Dell3, Peter Zeidler4, Olivier
Berne5, Carlo F. Manara6, Giulia Ballabio7, Jinyoung S. Kim8, John Bally9,
Javier R. Goicoechea10, Mari-Liis Aru6, Aashish Gupta6 and Anna Miotello6
1Astronomy Unit, School of Physics and Astronomy, Queen Mary University of
London, London E1 4NS, UK
2Department of Physics and Astronomy, Rice University, 6100 Main St - MS 108,
Houston, TX 77005, USA
3Department of Physics and Astronomy, Vanderbilt University, Nashville, TN
37235-1807, USA
4AURA for the European Space Agency (ESA), ESA Office, Space Telescope Science
Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA
5 Institut de Recherche en Astrophysique et Planétologie, Université de
Toulouse, CNRS, CNES, UPS, 9 Av. du colonel Roche, 31028 Toulouse Cedex 04,
France
6 European Southern Observatory, Karl-Schwarzschild-Strasse 2, 85748 Garching
bei München, Germany
7 Astrophysics Group, Imperial College London, Blackett Laboratory, Prince
Consort Road, London SW7 2AZ, UK
8 Steward Observatory, University of Arizona, 933 N. Cherry Ave, Tucson, AZ
85721-0065, USA
9 Center for Astrophysics and Space Astronomy, Astrophysical and Planetary
Sciences Department, University of Colorado, UCB 389 Boulder, Colorado 80309,
USA
10 Instituto de Física Fundamental (CSIC). Calle Serrano 121-123, 28006,
Madrid, Spain.
<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We present VLT/MUSE Narrow Field Mode (NFM) observations of a pair of disc-
bearing young stellar objects towards the Orion Bar: 203-504 and 203-506. Both
of these discs are subject to external photoevaporation, where winds are
launched from their outer regions due to environmental irradiation.
Intriguingly, despite having projected separation from one another of only
1.65″(660 au at 400 pc), 203-504 has a classic teardrop shaped “proplyd”
morphology pointing towards $\theta^{2}$Ori A (indicating irradiation by the
EUV of that star, rather than $\rm\theta^{1}$ Ori C) but 203-506 has no
ionisation front, indicating it is not irradiated by stellar EUV at all.
However, 203-506 does show [C i] 8727 Å and [O i] 6300 Å in emission,
indicating irradiation by stellar FUV. This explicitly demonstrates the
importance of FUV irradiation in driving mass loss from discs. We conclude
that shielding of 203-506 from EUV is most likely due to its position on the
observers side of an ionized layer lying in the foreground of the Huygens
Region. We demonstrate that the outflow HH 519, previously thought to be
emanating from 203-504 is actually an irradiated cloud edge and identify a new
compact outflow from that object approximately along our line of sight with a
velocity $\sim 130$ km s-1.

###### keywords:

stars: formation – ISM: H ii regions – ISM: kinematics and dynamics – stars:
jets – stars: protostars – planets and satellites: formation

††pubyear: 2023††pagerange: The VLT MUSE NFM view of outflows and externally
photoevaporating discs near the Orion Bar††thanks: Based on observations
collected at the European Southern Observatory under ESO programme
110.259E.001.–The VLT MUSE NFM view of outflows and externally
photoevaporating discs near the Orion Bar††thanks: Based on observations
collected at the European Southern Observatory under ESO programme
110.259E.001.

##  1 Introduction

With ever increasing numbers of exoplanets being discovered (Gaudi et al.,
2021), there is a need to understand how the planet formation process
contributes to exoplanet diversity. There is now overwhelming evidence that
planets form from circumstellar “planet-forming” discs, to the extent that we
can even directly observe young planets within them in some cases (Pinte et
al., 2018; Teague et al., 2018; Keppler et al., 2018; Haffert et al., 2019;
Pinte et al., 2022). Our most detailed studies of these have typically been
targeting large radius discs in the most nearby star forming regions, which
are the easiest to detect and resolve features in. Those have provided great
insight into the gas and dust composition and conditions, revealing
substructures that may possibly be associated with the planet formation
process, or planets themselves (e.g. ALMA Partnership et al., 2015; Dipierro
et al., 2015; Andrews et al., 2018; Zhang et al., 2018; Pinte et al., 2020;
Öberg et al., 2021). While those observations of nearby bright discs are
extremely important, recently there is a shift towards trying to study more
typical planet-forming discs. For nearby regions, this entails studying in
detail those of more typical radii($R\leq$50 au, Zhang et al., 2023; Miotello
et al., 2021, Miley et al. submitted, and the DECO ALMA large program:
2022.1.00875.L, PI: Cleeves). However, for a truly general understanding of
planet-forming disc evolution and planet formation we need to look to more
distant, more massive, stellar clusters.

The Sun is currently located in the middle of a supernova driven bubble that
probably induced our nearest ($\sim 140$ pc distant) instances of star
formation (for example, in Taurus and Lupus Zucker et al., 2022). These are
all low mass star forming regions, containing of order hundreds of stars, all
of which are $<10$ M⊙. The observed cluster mass function of $dN/dM\propto
M^{-2}$ implies that most stars form in much larger massive clusters (e.g.
Lada & Lada, 2003; Chandar, 2009; Fall et al., 2010). Massive star clusters
form from the clumpy/hierarchical collapse (or merger) of turbulent giant
molecular clouds (Lada & Lada, 2003; Portegies Zwart et al., 2010; Fukui et
al., 2021). Once massive stars form their radiative (photoionisaton/radiation
pressure) and mechanical (winds/supernovae) feedback heats and disperses the
cloud (Gritschneder et al., 2009, 2010; Dale & Bonnell, 2011; Dale et al.,
2012; Dale et al., 2014; Walch et al., 2012; Grudić et al., 2021; Grudić et
al., 2022; Guszejnov et al., 2022). The high stellar number density, dispersal
of the molecular cloud and high UV radiation field strengths all mean that the
star forming environment in which the planet-forming discs reside is very
different to lower mass star forming regions.

Massive stellar clusters have high stellar number densities, so the
probability of gravitational encounters between star-disc systems is more
likely (Pfalzner et al., 2005; Pfalzner, 2013; Pfalzner et al., 2018; Cuello
et al., 2019, 2020, 2023). Additionally, OB stars dominate the production of
UV photons by many orders of magnitude (for example a 40 M⊙ star emits around
ten trillion times more ionising photons than a Solar type star). One impact
of that is molecular cloud dispersal, which may influence the timescale that
discs can be fed by accretion streamers (e.g. Tang et al., 2012; Kratter &
Lodato, 2016; Yen et al., 2019; Pineda et al., 2020; Ginski et al., 2021;
Garufi et al., 2022; Valdivia-Mena et al., 2022; Huang et al., 2022; Gupta et
al., 2023). The much more established impact of the UV from massive stars is
that it leads to “external” photoevaporative winds being launched from the
outer parts of irradiated discs (for recent reviews see Winter & Haworth,
2022; Reiter & Parker, 2022). External photoevaporative winds are spatially
resolvable in relatively nearby sources ($D\sim 400$ pc), with teardrop shape
ionisation front morphologies now referred to as “proplyds” that point towards
the UV source responsible for irradiating them (e.g. O’Dell et al., 1993;
O’Dell & Wen, 1994; Henney & O’Dell, 1999; Bally et al., 2000; Ricci et al.,
2008; Miotello et al., 2012; Fang et al., 2016; Kim et al., 2016; Haworth et
al., 2021; Winter & Haworth, 2022). Proplyd mass loss rates are empirically
inferred to be sufficiently high that it would result in rapid depletion and
truncation (Henney & O’Dell, 1999; Winter et al., 2019; Winter & Haworth,
2022) and models support these high mass loss rates (Richling & Yorke, 2000;
Adams et al., 2004; Facchini et al., 2016; Haworth et al., 2018; Haworth &
Clarke, 2019). Recent models predict that external photoevaporation could
impact planet formation anywhere in the disc, by reducing the reservoir of
solids that can radially drift and resupply the inner disc with material (Qiao
et al., 2023), as well as affect the migration and growth of giant planets.
The discs of proplyds have also been subject to detailed study to understand
the impact of the external winds upon them, finding that they are typically
more compact, less massive and even chemically distinct from discs in low mass
star forming regions (e.g. Mann et al., 2014; Eisner et al., 2018; Boyden &
Eisner, 2020, 2023).

Because the impact of external photoevaporation on discs is expected to be so
significant, it makes the connection to the wider star formation process
particularly important. Stars form embedded and hence shielded from external
photoevaporation until they dynamically emerge from the molecular cloud, or
the cloud itself is dispersed by feedback. This has been shown to be important
for understanding disc properties (masses, radii, lifetimes) in clusters (Qiao
et al., 2022; Wilhelm et al., 2023; Concha-Ramírez et al., 2023), and the
shielding timescale is predicted to have a significant impact on the resulting
planetary architecture (Qiao et al., 2023).

In this paper we study two externally irradiated discs coincident with the
Orion Bar. One of these is a classic example of an externally photoevaporating
disc, with an extreme ultraviolet (EUV) excited teardrop shaped ionisation
front. The other is somehow shielded from EUV irradiation, has a wind driven
purely by far-ultraviolet (FUV). This pair of objects is valuable for
understanding the role of external photoevaporation outside of the area of
most intense feedback in the Orion Nebular Cluster, and the role that
intervening material plays in changing the UV field incident upon discs.

##  2 The Orion Nebula and the Orion Bar

At a distance of $\sim 400\,$pc (Kounkel et al., 2017, 2018), the Orion Nebula
Cluster is the nearest and best studied massive star forming region. The
primary source of UV photons in the region is from the Trapezium stars, and of
those most notably the O6V star $\theta^{1}$ Ori C (O’Dell et al., 2017b).
Figure 1 shows a Hubble Space Telescope 3-colour ([N ii], H $\alpha$, [O iii])
image of the overall region, with the main features labeled, including the
location of the main field of view (FOV) that we consider in this paper. This
HST image is made from the mosaic of emission-line filters (red-[N II] 6583 Å,
green-H$\,\alpha$, blue-[O III] 5007 Å) first presented in O’dell & Wong
(1996), but reprocessed and improved by Jeremy R. Walsh and Richard Hooks of
the European Southern Observatory.

Figure 1: A 3-colour HST [N ii] (red) H $\alpha$ (green) and [O iii] 5007 Å
image image of the Trapezium (upper right) and Orion Bar (also known as the
Orion Bright Bar). The FOV that we study with VLT MUSE is the boxed region in
the south western Orion Bar (which O’Dell et al. (2017a) call the SW Bright
Bar). The orientation of ionisation fronts associated with the most prominent
proplyds is denoted by the white arrows. The main UV source in the region is
$\theta^{1}$ Ori C in the Trapezium, but $\theta^{2}$ Ori A is also clearly
irradiating discs and has been shown to dominate UV irradiation of parts of
the cloud south east of the Orion Bar O’Dell et al. (2017b). Note that North
is up in this image.

The overall ionisation structure around $\theta^{1}$ Ori C is as follows. The
innermost photoionised region contains ionised hydrogen (ionisation energy
13.6 eV) and singly ionised helium (ionisation energy 24.59 eV). There is no
doubly ionised helium zone because $\theta^{1}$ Ori C does not emit enough
sufficiently energetic ($>54.4$ eV) photons. This singly ionised helium zone
is most easily traced by the [O iii] 5007 Å emission line from doubly ionised
oxygen. Moving outwards through the photoionised gas, all of the photons
capable of ionising helium are absorbed, resulting in a zone of neutral helium
but ionised hydrogen that is most easily traced by [N ii] 6583 Å emission. At
the distance where all of the hydrogen ionising photons are used in
maintaining photoionisation balance in the H ii region is a relatively sharp
ionisation front, beyond which is a photodissociation region (PDR).

The brightest area of the Orion Nebula is historically called the Huygens
Region. Quantitative analysis of emission lines there has led to the
interpretation of the Huygens Region being a concave surface with
irregularities(Wen & O’Dell, 1995; Henney et al., 2005; O’Dell et al., 2009,
2017b).

In this picture, the Trapezium stars are near the more distant side of the
concave surface, which slopes closer to the observer towards a bright, linear
southeast feature. The main ionisation front (MIF) associated with this
feature is where the hydrogen ionisation front is along the line of sight. In
the foreground is a “nearer ionised layer” (NIL O’Dell et al., 2020). The
density and temperature structure of the Huygens Region has been well studied
(e.g. Weilbacher et al., 2015; McLeod et al., 2016; O’Dell et al., 2017a).

The Huygens region of the Orion Nebula contains multiple elongated bar
structures (O’Dell et al., 2017a) but the most famous, around
$112^{\prime\prime}$ to the south east of the Trapezium is the NE-SW tending
Orion Bar (alternatively referred to as “the Bright Bar" in many visual
studies, see Figure 1). In this report we will subsequently call it the “Bar”.
The Bar is a feature caused by the MIF of the central Orion Nebula being
tilted so much that it is viewed edge-on. This allows a stacked view of the
structure in the underlying PDR. The Bar has been the subject of a number of
1D modelling efforts (e.g. Bregman et al., 1989; Tauber et al., 1994; Jansen
et al., 1995; Walmsley et al., 2000; Cuadrado et al., 2015), as well as models
and observations exploring the importance of clumpiness and 3D structure in
PDRs (Andree-Labsch et al., 2017). Goicoechea et al. (2016) find that
considering the Orion bar as a static slab may be insufficient, since it is
being compressed and ablated by the high pressure photoionised gas and/or
photoevaporative compression by the rocket effect. Habart et al. (2023) also
demonstrate the complicated sub-structured nature of the Orion Bar.

The region that we study in this paper contains two externally
photoevaporating young stellar objects (YSOs), 203-506 and 203-504, at the
boundary of the Bar. Their coordinates are given in Table 1. Both objects have
an associated outflow: HH 519 and HH 520, respectively. HST has already
demonstrated that 203-504 is a proplyd that is directed towards $\theta^{2}$
Ori A (Bally et al., 2000; O’Dell et al., 2017b). 203-506 was included in the
Herschel study of Champion et al. (2017), Habart et al. (2023) recently
detected vibrationally excited (v=1-0 S(1) at 2.12 $\mu$m) H2 emission, and
Berne et al. (submitted) have most recently used JWST (Berné et al., 2022) and
ALMA to study 203-506.

Our objective has been to use the European Southern Observatory (ESO) integral
field spectrograph MUSE on the Very Large Telescope (VLT), which in the Narrow
Field Mode (NFM) has spatial resolution comparable to HST and spectral
coverage over the entire optical range to understand more about these
apparently very different externally photoevaporating discs, which are seen in
projection close to the boundary of the Huygens Region, just beyond the peak
brightness of the Bar. The multiplicity of emission-lines in the visual region
were expected to give a better depiction of the environment of both 203-504
and 203-506 and a refinement of the physical conditions within 203-506.

System | RA | DEC |
---|---|---|---
203-504 (HH 519) | 05 35 20.26 | -05 25 04.05 |
203-506 (HH 520) | 05 35 20.32 | -05 25 05.55 |
Table 1: The coordinates of the YSOs 203-504 and 203-506 (Ricci et al., 2008).

###  2.1 The distance to $\theta^{2}$ Ori A

Although not the dominant UV source in the Orion Nebular Cluster (ONC), the
O9.5V star $\theta^{2}$ Ori A plays an important role in this paper. It is
located to the south east of the Bar (it is labeled in Figure 1). The details
of our interpretation are sensitive to the distance of $\theta^{2}$ Ori A
relative to the rest of the cluster. The distance to the ONC is generally
inferred to be around 390-400 pc (Jeffries, 2007; Kounkel et al., 2017, 2018;
Großschedl et al., 2018). The Gaia DR3 and DR2 distances to $\theta^{2}$ Ori A
are 336 and 450 pc respectively (Gaia Collaboration et al., 2016, 2018, 2022)
so are wildly inconsistent with one another, suggesting that the current Gaia
estimates are affected by some factor such as its multiplicity and/or larger
uncertainties in the parallax of bright sources (Schulz et al., 2006) or
nebulosity of the region. Given its importance and our inability to use a Gaia
distance in this case we review the other body of evidence on its location,
which supports $\theta^{2}$ Ori A being at a similar distance as the main
cluster.

The strongest argument lies in the two-pronged study of O’Dell et al. (2017b)
which addresses the sources of ionization of the Orion Nebula and its MIF.
There they established that the bright ionized rims of the proplyds near
$\rm\theta^{1}$ Ori C were oriented towards that star. However, when
proceeding to the SE, the proplyd bright rims transition to orientation
towards both $\rm\theta^{1}$ Ori C and $\rm\theta^{2}$ Ori A, until becoming
exclusively oriented towards $\rm\theta^{2}$ Ori A (see Figure 1, where white
arrows denote the proplyd orientations). They also examined the line ratios of
the brightest lines in both the visual and infrared windows, using the cloudy
code (Ferland et al., 1998; Ferland et al., 2013, 2017) and the known
radiation properties of the stars. They found the MIF in the Huygens Region to
be photoionized by $\rm\theta^{1}$ Ori C, but to the southeast of the Bar,
$\rm\theta^{2}$ Ori A was the source of ionization.

The radial velocity of $\rm\theta^{2}$ Ori A is compatible with its membership
in the ONC. Although the SIMBAD data-base gives this as 19 km s-1 citing an
early low resolution study (Henroteau & Henderson, 1920), examination of that
paper shows that although the star was recognized as a spectroscopic binary,
the results were confused by the presence of emission lines in the spectra. A
definitive value for the system velocity was determined from the orbital
solution of Abt et al. (1991) to be 24.8$\pm$2.1 km s-1. This is essentially
the same as the most frequent velocity of the ONC stars of 25$\pm$2 km s-1 in
the study of Sicilia-Aguilar et al. (2005). This agreement would not be
expected if $\rm\theta^{2}$ Ori A was not a cluster member. Note that all
velocities in this paper are heliocentric, for LSR subtract 18.2 km s-1, which
is arrived at under the usual assumption that the Sun is moving at 20 km/s
towards 18:00:00 +30:00:00 (epoch 1900).

Overall the above arguments all point towards $\theta^{2}$ Ori A not being
significantly in the foreground, rather, being close enough to the nebula to
play a role in its photoionization processes.

##  3 Observations and data reduction

We were awarded 1 hour of Directors Discretionary Time (DDT) on VLT MUSE (PI:
Haworth, run ID:110.259E.001). MUSE (Multi Unit Spectroscopic Explorer) is the
visual wavelength integral-field spectrograph on the VLT mounted on UT4 (Bacon
et al., 2010). Our observations were a single MUSE narrow field mode (NFM)
pointing with adaptive optics (AO). The observations were taken in service
mode on the 20th November 2022 with 4 integrations of 11 minutes, giving 44
minutes on source. 4 raw datasets were provided in the observing block at 90
degree orientations to one another. These were reduced separately and then
combined to make a single datacube. We reduced the data using the MUSEpack
(Zeidler et al., 2019; Zeidler, 2019) python wrapper for the standard esorex
(v2.8.7) framework (Freudling et al., 2013). This utilised the reduced master
calibration files provided by ESO.

The FOV of the data is approximately $7.42\arcsec$ per side with an effective
angular resolution of $0.0921\arcsec$, corresponding to spatial scales of 2968
au and 36.8 au respectively. The wavelength range spans 4750-9350 Å and the
spectral resolving power over that range varies from $R=2000-4000$,
corresponding to velocity resolution of $\sim 75-150$ km s-1. We hence do not
generally expect to resolve the kinematics of our targets, except possibly in
the case of high velocity outflows aligned along the line of sight.

The sky subtraction in the reduction of MUSE data uses either a dedicated sky
observation, or the darkest part of the image. However, for bright and
extended targets such as the Orion star forming region it is not possible to
obtain an appropriate sky measurement because the telescope has to move by
many degrees to find an appropriately dark field. In regions with high
background, such as an H  ii region, using the darkest part of an image may
lead to over subtraction of key (nebular) lines of interest. To resolve this,
we use the modified sky subtraction approach introduced by Zeidler et al.
(2019), which is included in MUSEpack. This assumes that the H  ii region
lines and continuum dominate over the sky, and uses OH and O2 to determine and
subtract the sky contribution. Even though this method may lead to a small
overestimation of the source line fluxes (originated in line emission of
Earth’s atmosphere), it has been shown that overall this method performs
better in terms of line profiles, flux ratios, and absolute wavelength
(velocity) measurements (Zeidler et al., 2019).

To study lines we extract data $\pm 3$Å from the nominal line centre. The
continuum subtraction uses nearby line free channels, typically 10 Å from the
line centre of interest so long as it is clear of other line emission.

To correct for reddening we follow O’Dell & Harris (2010); Weilbacher et al.
(2015); Reiter et al. (2019) and use the Balmer decrement to estimate the line
of sight extinction. That is, we assume an intrinsic H$\,\alpha$/H $\beta$
ratio of 2.89 (Storey & Hummer, 1995; Osterbrock & Ferland, 2006) and use the
pyneb (Luridiana et al., 2015) reddening correction routine with a Cardelli et
al. (1989) and Blagrave et al. (2007) extinction law and assuming an optical
total-to-selective extinction ratio $R_{V}=5.5$ (Cardelli et al., 1989;
Baldwin et al., 1991; Osterbrock et al., 1992; Greve et al., 1994; Blagrave et
al., 2007; Weilbacher et al., 2015; Fang et al., 2021). That extinction law
and higher than standard value of $R_{V}$ for the diffuse ISM was found to be
more appropriate for the ONC in a VLT MUSE study of hundreds of stars in the
ONC by (Fang et al., 2021). A value of $R_{V}=5.5$ was also adopted in the
MUSE study of the ONC by Weilbacher et al. (2015).

##  4 Results

Figure 2: A 3 colour composite image from our dataset consisting of [C i] 8727
Å (red), H$\,\alpha$ (green) and [O i] 8446 Å (blue). The two disc-bearing
YSO’s are 203-504 and 203-506. To the upper left in green is the H ii region
that is photoionised by the Trapezium stars (predominantly the O6.5V star
$\theta^{1}$ Ori C). The O9.5V star $\theta^{2}$ Ori A is towards the left of
the FOV. Blue traces Lyman $\beta$ irradiated gas, which is predominantly at
the interface with the H ii region, but also traces an irradiated cloud near
the line of sight location of the two YSOs, which based on its shape we call
the Delta ($\Delta$) Feature.

###  4.1 Overview

We begin with an overview of our dataset and its key components. Figures 2 and
3 show colour composites, and Figure 4 shows a gallery of individual continuum
subtracted emission lines.

Figure 2 combines [C i] 8727 Å (red), H $\alpha$ (green) and [O i] 8446 Å
(blue). We will later argue in Appendix B that the [C i] 8727 Å line traces
recombinations in FUV irradiated gas, and in our field is predominantly
tracing the circumstellar material of 203-506. H $\alpha$ is a well
established tracer of photoionised gas. The [O i] 8446 Å traces Lyman $\beta$
irradiated gas because H i Lyman $\beta$ radiation ($\lambda=1025.72$ Å) has a
very similar wavelength to the 2$\rm p^{3}$ 3d $\rm{}^{3}D^{o}_{3}$ ground
state of atomic oxygen ($\lambda=1025.76$ Å). The [O i] 8446 Å line is a
forbidden line emitted as part of the decay following resonant absorption of
the Lyman $\beta$ as described by Bowen (1928) and also detailed in Figure 4.7
of Osterbrock & Ferland (2006). Usually the 8446 Å line is weak in Huygens
Region spectra, but it is often strong in our spectra. We have adopted the
Bowen (1928) process described above as the source of the emission, with its
requirement that Oo is being irradiated by hydrogen Ly$\beta$ radiation (see
also Osterbrock & Ferland, 2006).

The most striking features in Figure 2 are the two externally irradiated
protoplanetary discs, 203-504 and 203-506. 203-504 exhibits the classic
teardrop shaped H $\alpha$ morphology of proplyds, pointing towards
$\theta^{2}$ Ori A. 203-506 is distinct, showing no clear ionisation front. We
will discuss the very different external photoevaporation of these discs
further in section 4.2.

Figure 3 again shows the [C i] 8727 Å (red) and [O i] 8446 Å (blue), but the
green in this case is the [Fe ii] 8617 Å line. [Fe ii] emission is well-
established as a tracer of shock-excited gas in protostellar jets (e.g.
Pesenti et al., 2003; Hartigan et al., 2004; Agra-Amboage et al., 2011;
Giannini et al., 2013) and has also been observed in a number of externally-
irradiated jets in H ii regions (e.g. Reiter & Smith, 2013; Reiter et al.,
2015; Reiter et al., 2016). In particular, Kirwan et al. (2023) used several
[Fe ii] lines observed with MUSE NFM to trace the properties of the jet from
the 244-440 proplyd, also irradiated from $\theta^{2}$ Ori A. In our dataset
[Fe ii] traces outflows from both 203-506 (HH 520) and 203-504 (HH 519), which
we discuss in more detail in 4.3.

To the upper right in Figure 2 is the Huygens Region, photoionised by the
Trapezium stars and bright in H $\alpha$.The sharp boundary indicates where
the MIF is nearly along the line-of-sight. The MIF then flattens into the
lower surface-brightness region we call the Bar Shoulder and continues beyond
to the southeast with progressively lower surface-brightness. The low
ionization boundary of the MIF is visible in Lyman $\beta$ irradiated [O i]
8446 Å. To the south-east, beyond the Bar ionisation front is the Delta
Feature, that is enhanced in [O i] 8446 Å. 203-506 is along the line of sight
of the Delta Feature, appearing in silhouette in many emission-lines (see
Figure 4) and 203-504 is located slightly to the northwest of the northern
vertex of the Delta Feature. We discuss the possible role that this feature
plays in the context of the externally photoevaporating discs in section 5.2.
Table 2 provides a non-exhaustive summary of bright lines in our dataset and
the features that they trace.

Finally, within the H ii region in the upper right of Figure 2 are three small
dark clumps, which are instrumental artefacts. We discuss these in Appendix A.

Figure 3: A 3 colour composite image from our dataset consisting of [C i] 8727 Å (red), [Fe ii] 8617 Å (green) and [O i] 8446 Å (blue). The [Fe ii] traces the bipolar outflow, HH 520, from 203-506 and the compact outflow HH 519 from 203-504.  Figure 4: A gallery of some key continuum-subtracted lines towards 203-504/203-506 in the Bar. The upper left panel is [C i], which traces FUV irradiated material in the discs. The upper right panel is [Fe ii], which traces shock excited gas from jets/outflows. Centre-left and centre-right are H$\,\alpha$ and [S ii] which trace EUV-irradiated gas. The lower panels are [O i] 6300 Å (left) which can trace collisionally and shock excited gas as well as regions of OH dissociation, and the [O i] 8446 Å line (right) which traces Lyman $\beta$ irradiated gas. Line | Wavelength | Primarily Traces | Notes
---|---|---|---
| (Å) |  |
H $\beta$ | 4861 | The H  ii region. Proplyd 203-504. The irradiated Delta Feature. |
$\left[\textrm{O\,{\,{ii}i}}\right]$ | 5007 | The H  ii region. |
$\left[\textrm{N\,{ii}}\right]$ | 5755 | The H  ii region. Proplyd 203-504. The irradiated Delta Feature. |
$\left[\textrm{OI}\right]$ | 6300 | HH 519. 203-506 disc and inner part of HH 520. Delta Feature and Bar PDR. |
$\left[\textrm{S\,{\,{ii}i}}\right]$ | 6312 | The H  ii region. Proplyd 203-504 (faint). |
$\left[\textrm{O\,{i}}\right]$ | 6363 | HH 519. 203-506 disc and inner part of HH 520. Delta Feature and Orion Bar PDR. |
H $\alpha$ | 6582.8 | The H  ii region. Proplyd 203-504. The irradiated Delta Feature. | 203-506 in silhouette
He i | 6678 | The H  ii region. Proplyd 203-504 | 203-506 in silhouette
$\left[\textrm{N\,{ii}}\right]$ | 6548 | The H  ii region. Proplyd 203-504. The irradiated Delta Feature. |
$\left[\textrm{N\,{ii}}\right]$ | 6583 | The H  ii region. Proplyd 203-504. The irradiated Delta Feature. |
$\left[\textrm{S\,{ii}}\right]$ | 6716 | The H  ii region. Part of HH 520 | 203-506 in silhouette
$\left[\textrm{S\,{ii}}\right]$ | 6732 | The H  ii region |
He i | 7065 | The H  ii region. Proplyd 203-504 | 203-506 in silhouette
$\left[\textrm{Fe \,{ii}}\right]$ | 7155 | The outflows, HH 519 and HH 520 |
$\left[\textrm{Fe \,{ii}}\right]$ | 7172 | The outflows, HH 519 and HH 520 |
$\left[\textrm{O\,{i}}\right]$ | 7255 | 203-504 proplyd. Delta Feature and Bar PDR. | 203-506 in silhouette
He i | 7282 | The H  ii region. Proplyd 203-504 | 203-506 in silhouette
$\left[\textrm{Ca \,{ii}}\right]$ | 7291 | The outflows, HH 519 and HH 520 |
$\left[\textrm{Ar \,{\,{ii}i}}\right]$ | 7751 | The H  ii region. Proplyd 203-504 | 203-506 in silhouette
$\left[\textrm{O\,{i}}\right]$ | 8446 | 203-504 proplyd. Delta Feature and Bar PDR. | 203-506 in silhouette
$\left[\textrm{Fe \,{ii}}\right]$ | 8617 | The outflows, HH 519 and HH 520 |
$\left[\textrm{C}\,\textsc{i}\right]$ | 8727 | Disc/wind of 203-504 and 203-506 |
$\left[\textrm{S\,{\,{ii}i}}\right]$ | 9069 | The H  ii region. Proplyd 203-504. |
Table 2: A non-exhaustive summary of lines detected in our MUSE dataset, and
the features they trace.

###  4.2 Two extremely different photoevaporating discs

Our FOV contains two protoplanetary discs subject to external
photoevaporation. Their projected on-sky separation is only around 660 au, but
the nature of their external photoevaporation is considerably different and
they are probably not located physically close to one another.

####  4.2.1 203-504

203-504 has a classic teardrop shaped proplyd morphology that is clear in
ionisation front tracers such as H $\alpha$. The 203-504 proplyd points
towards $\theta^{2}$ Ori A, implying that is the UV source predominantly
responsible for driving the mass loss (as already pointed out by O’Dell et
al., 2017b). In addition, there is [C i] 8727 Å and [O i] 8446 Å emission from
FUV/Lyman $\beta$ irradiated gas (see Figure 4). 203-504 is hence being
exposed to both FUV and EUV radiation.

In ionisation equilibrium, the extent of the ionisation front in a proplyd
depends on the ionising flux incident upon the proplyd and the integrated
recombination rate from the UV source down to the I-front. Assuming that the
flow diverges spherically from the I-front and that a spherical flow from the
disc into the I-front has to satisfy mass conservation, one can estimate the
required mass loss rate to place the I-front at its observed radius. For
example, for a given incident ionising flux, decreasing the mass-loss rate
would make the ionisation front radius smaller (e.g. Johnstone et al., 1998;
Clarke & Owen, 2015). Specifically, for a proplyd at a projected separation
$d$ from a UV source producing ionising photons at a rate $\dot{N}_{ly}$ per
second and with mass-loss rate $\dot{M}_{\textrm{wind}}$, the ionisation front
radius is

$R_{\mathrm{IF}}\approx
1200\,\left(\frac{\dot{M}_{\mathrm{wind}}}{10^{-8}\,M_{\odot}\,\mathrm{yr}^{-1}}\right)^{2/3}\\\
\times\left(\frac{\dot{N}_{\mathrm{Ly}}}{10^{45}\,\mathrm{s}^{-1}}\right)^{-1/3}\left(\frac{d}{1\,\mathrm{pc}}\right)^{2/3}\,\mathrm{au}$
(1)

(Winter & Haworth, 2022).

The projected separation from 203-504 to $\theta^{2}$ Ori A is $d=0.077$ pc
assuming a distance of 400 pc. For the ionising luminosity of $\theta^{2}$ Ori
A we use the value from Table 2.3 of Osterbrock & Ferland (2006) of
$3.63\times 10^{48}$ photons s-1. We measure an ionisation front radius for
203-504 of 139 au, which combined with the values above corresponds to a mass
loss rate of $\sim 3.1\times 10^{-7}$ M⊙ yr-1. In making this estimate we have
assumed based on the I-front morphology that only $\theta^{2}$ Ori A is
driving the mass loss. We will discuss the possible relative locations of the
evaporating discs/UV sources and the ISM in section 5.2.

As discussed, the EUV field sets the ionisation front radius, but the FUV
radiation is expected to be what actually drives the wind launching and sets
the mass-loss rate (Johnstone et al., 1998; Adams et al., 2004; Haworth et
al., 2018, Haworth et al. submitted). Assuming that the FUV incident upon
203-504 is also dominated by $\theta^{2}$ Ori A and that $\theta^{2}$ Ori A
and 203-504 are the same distance from the observer, integrating over a Kurucz
(1979) model spectrum suggests an FUV field incident upon 203-504 of
approximately $8\times 10^{4}$G0. Other proplyds at a comparable FUV field
strength are seen for in NGC 2024, which are irradiated by an O8V star and
also have mass loss rates of order $10^{-7}$ M⊙ yr-1 (Haworth et al., 2021).
The majority of proplyds in the ONC are in the vicinity of $\theta^{1}$ Ori C
irradiated by $\geq 10^{6}$ G0 UV fields with correspondingly higher mass loss
rates (regularly $\geq 10^{-6}$ M⊙ yr-1).

####  4.2.2 203-506

The other young stellar object, 203-506 has a drastically different
morphology, not seen in any other proplyds. The main feature that is unusual
in this system is that it has no ionisation front, evidence that it is not
being irradiated by EUV photons. It is however being irradiated by the FUV,
which is responsible for driving the [C i] 8727 Å (see Appendix B) and the [O
i] 6300 Å line by OH dissociation (see section 5.4).

203-506 is bright in [O i] 6300 Å and [C i] 8727 Å, but appears in silhouette
in H$\alpha$, [O iii] 5007 Å, [N ii] 6583 Å, [O i] 8446 Å, and [S ii] 6732 Å.
Reference to images made immediately after installation of the WFPC2 on the
HST show the same properties in the superior H$\alpha$ F656N and [N ii] F658N
images. The fact that the silhouettes are similar in all ions (except when
203-506 is in emission) means that 203-506 lies well into the foreground of
the MIF.

###  4.3 Jets/outflows

Our FOV contains two Herbig-Haro (HH) objects: HH 519 and HH 520 associated
with 203-504 and 203-506 respectively. They were originally identified by
Bally et al. (2000) as a pair of oppositely directed monopolar outflows in
narrowband [O i] 6300 Å images, as well as H $\alpha$ in the case of HH 519.

In our data, HH 520 (emanating from 203-506) is clearly observable in jet
tracing (Reiter & Smith, 2013; Reiter et al., 2015; Reiter et al., 2016;
Kirwan et al., 2023) [Fe II] emission, specifically the 8617 Å a4F9/2-a4P5/2
line. Of note is that we also detect fainter, oppositely directed, [Fe ii]
emission from HH 520 (see Figure 3 and the upper right panel of Figure 4).
MUSE reveals the eastern limb of the jet for the first time, demonstrating
that this is a bipolar flow, but with the south eastern component a factor
5-10 fainter than the previously identified north western component.

Figure 5: A continuum-subtracted position-velocity diagram in a cut across the
YSO 203-504, with the reference velocity corresponding to the [Fe ii] 8617 Å
line. The velocity of the outflow HH 519 is resolved at $\sim 130\,$km s-1.
The lower panel shows the location of the cut (compare with the upper right
panel of Figure 5).

Our data for HH 519 (for example in Figure 4) demonstrates that the H $\alpha$
previously identified as the outflow actually traces the northeast bright rim
of the Delta Feature. The underlying cloud is particularly clear in the [O i]
6300 Å and [O i] 8446 Å emission in the lower panels of Figure 4. 203-504 is
coincidentally located close to the northern vertex of that irradiated Delta
Feature along the line of sight, making it appear as though the elongated
bright rim originated from it. However, the high angular resolution of MUSE
NFM (0.0921”) also demonstrates that 203-504 is completely disconnected from
the irradiated Delta Feature (particularly in [O i] 8446Å ).

Our data also reveal a compact and kinematically resolved high velocity
outflow feature in [Fe ii] 8617Å associated with 203-504. Because the outflow
is associated with the same YSO the name “HH 519” is retained for this new
outflow associated with 203-504 (it is labeled in Figure 3 and in the upper
right panel of Figure 4). In Figure 5 we make a cut (a single pixel wide)
across the HH 519 jet axis in continuum subtracted [Fe ii] emission and show
the position-velocity diagram. The outflow appears to be bipolar, with each
outflow limb extending $\sim 0.4^{\prime\prime}$ and line of sight velocity
$\sim 130\,\textrm{km}\,\textrm{s}^{-1}$ that is consistent with velocities
measured in other HH jets (e.g. Hartigan et al., 2001, 2005; Hartigan & Morse,
2007; Reiter et al., 2017, 2022; Kirwan et al., 2023). The compact size and
fact that it is spectrally resolved suggests that HH 519 is tilted toward the
line of sight by $\gtrsim 45^{\circ}$.

##  5 Discussion

###  5.1 The Irradiation of our FOV

####  5.1.1 Strong Line ratios argue for irradiation by $\rm\theta^{1}$ Ori C
and/or $\rm\theta^{2}$ Ori A

O’Dell et al. (2017b) used line ratios to determine the source of UV
irradiation at the MIF across the Huygens region, and were able to distinguish
when $\theta^{1}$ Ori C and $\theta^{2}$ Ori A were the dominant source (with
one field that they studied likely receiving significant contributions from
both sources). To compare with O’Dell et al. (2017b), three small sub-samples
were taken. Within the Delta Feature the sub-sample was
0$\aas@@fstack{\prime\prime}$52$\times$0$\aas@@fstack{\prime\prime}$52 centred
0$\aas@@fstack{\prime\prime}$78 east and 0$\aas@@fstack{\prime\prime}$15 south
of 203-506. Outside of the Delta Feature an east sub-sample was
0$\aas@@fstack{\prime\prime}$69$\times$1$\aas@@fstack{\prime\prime}$25 centred
2$\aas@@fstack{\prime\prime}$38 east and 0$\aas@@fstack{\prime\prime}$37 north
of 203-506 and a west sub-sample of
0$\aas@@fstack{\prime\prime}$80$\times$0$\aas@@fstack{\prime\prime}$61 centred
1$\aas@@fstack{\prime\prime}$39 west and 0$\aas@@fstack{\prime\prime}$07 north
of 203-506. The averages of these three sub-samples are shown in Table 3, and
compared also with the other measurements in the Huygens Region. The values in
our field lie between those in the central Huygens Region and those to the SE
of the Bar, as one would expect given our field lies between the two. O’Dell
et al. (2017b) determined that their sample E (OKF-E) is likely subject to
significant irradiation by both $\theta^{2}$ Ori A and $\theta^{1}$ Ori C. Our
FOV is most similar to OKF-E, suggesting both stars may also play a role here,
though with perhaps a fractionally more significant contribution from
$\theta^{1}$ Ori C.

Table 3: Huygens Region and Delta Line Ratios. OKF refers to a sample
presented in O’Dell et al. (2017b)

Sample | I([O iii] 5007)/I(H$\beta$) | I([N ii] 6583)/I(H$\alpha$)
---|---|---
Esteban et al. (2004)a | 3.84 | 0.13
OKF-Aa | 3.73 | 0.12
OKF-B a | 3.58 | 0.14
Ave 203-506b | 2.16$\pm$0.04 | 0.23$\pm$0.01
OKF-Ec | 1.96 | 0.18

aSamples near $\rm\theta^{1}$ Ori C.
bAverage of Delta Feature and two adjacent samples.
cSample in O’Dell et al. (2017b) nearest 203-504/6, which they determined was
affected by both $\theta^{2}$ Ori A and $\theta^{1}$ Ori C.

####  5.1.2 Comparison of line ratios over a wide range of ionization energies
argue for both $\rm\theta^{1}$ Ori C and $\rm\theta^{2}$ Ori A

We additionally compared the average line fluxes over a bigger range of
ionization in our samples described above to those in the vicinity of
$\theta^{1}$ Ori C in Table 4. The Inner Delta values are from the “Ave
203-506” samples in Table 3 and the Average Background samples are from
outside the Delta Feature. These fluxes are defined relative to
$F($H$\,\beta)=100$. We see that the high ionisation [O iii] 5007 Å line is
much weaker than in the vicinity of $\theta^{1}$ Ori C, while the intermediate
ionisation [N ii] 6853 Å line is stronger and the low ionisation [O i] 6300 Å,
[S ii] 6732 Å and [O i] 8446 Å lines are all much stronger. This is consistent
with the known pattern of decreasing ionization with increasing distance from
$\rm\theta^{1}$ Ori C O’Dell & Harris (2010). The average background and
fluxes in the inner part of the Delta Feature are similar, with the exception
of the [O i] 8446 Å line, which is elevated by over a factor 2 compared to the
background and by over a factor 3 compared to the vicinity of $\theta^{1}$ Ori
C, which as we discuss in more detail in appendix C.

Both this section and its predecessor both point towards the radiation field
in our data not being exceptional compared to other samples in the region,
i.e. that both stars can play a role in irradiating our FOV.

The presence of the distinct boundaries on the Delta Feature therefore point
to genuine structure, such as density enhancements and/or a raised clump in
the shoulder region of the Bar’s MIF.

Table 4: Line fluxes near $\theta^{1}$ Ori C, ”outside” and “inside” the Delta Feature. The fluxes are defined relative to the H$\,\beta$ line with a value of 100.  Ion,Line | Near $\rm\theta^{1}$ Ori Ca | Average Background | Delta Feature
---|---|---|---
H$\beta$ | 100 | 100 | 100
[O iii] ,5007 Å | 384 | 216 | 220
[O i] ,6300 Å | 0.71 | 1.0 | 1.6
H$\alpha$ | 287 | 293 | 293
[N ii] , 6583 Å | 38 | 70 | 65
[S ii] , 6732 Å | 3.5 | 7.8 | 8.6
[O i] , 8446 Å | 0.88 | 1.3 | 2.8
aFrom Esteban et al. (2004). |  |  |

####  5.1.3 Information from the Delta Feature

The Delta Feature has three edges, one in the north east, one in the north
west and one in the south. The NE and NW boundaries of the Delta Feature are
more sharply defined than the south boundary. They have the appearance of
ionization fronts viewed edge-on and are bright in [O i] 6300 Å, [S ii] 6732 Å
[N ii] 6583 Å, and [O i] 8446 Å. In H$\alpha$ and H$\beta$ they are more
diffuse, but are bright and in about the same positions. None of the edges of
the Delta Feature are associated with changes in [C i] 8727 Å or [O iii]
5007Å. These characteristics all point towards the rim of the Delta Feature
being irradiated by some EUV radiation with no significant radiation above
24.6 eV (due to the lack of the helium ionised zone tracing [O iii]). The NW
boundary is approximately perpendicular to a line towards $\theta^{2}$ Ori A
and the NW boundary is approximately perpendicular to a line towards
$\theta^{1}$ Ori C, suggesting that the NE and NW edges of the Delta Feature
are being irradiated by the two different UV sources. This is consistent from
the two preceding sections concluding that our FOV is likely to be influenced
by both $\rm\theta^{1}$ Ori C and $\rm\theta^{2}$ Ori A. The presence of the
distinct boundaries on the Delta Feature point to genuine structure, such as a
density enhancement and/or a raised clump in the shoulder region of the Bar’s
MIF.

###  5.2 The locations of the externally photoevaporating discs

Figure 6: A cartoon schematic (not to scale) of the Huygens region and the key
features in our field. $\theta^{1}$ Ori C sits within a concave surface. Where
that surface orients along the line of sight is the Bar ionisation front. Our
field, with the Delta Feature and the two evaporating discs, 203-504 and
203-506, lies in the direction of the shoulder region, where $\rm\theta^{1}$
Ori C’s radiation begins to be shadowed. The nearer ionised layer (NIL) is
somewhere in the foreground of $\theta^{1}$ C but the exact distance is not
well known, particularly near $\theta^{2}$ Ori A (O’Dell et al., 2020).

In section 4.2 we discussed that 203-504 and 203-506, with a projected
separation of 660 au on the sky, appear to be being subjected to external
photoevaporation by very different external UV radiation fields. 203-504 is
being irradiated by both direct EUV and FUV radiation, whereas 203-506 is
solely being irradiated by the FUV and is hence somehow shielded from the EUV
from $\theta^{2}$ Ori A and $\theta^{1}$ Ori C. Here we discuss the possible
structure of the region that could give rise to this scenario. To do this, we
can draw on existing information about the structure of the nebula from the
Huygens Region, past the Bar, and into the Extended Orion Nebula (EON).

####  5.2.1 Location of key stars and features that block portions of their
radiation

Figure 6 shows a cartoon schematic cross section from near $\rm\theta^{1}$ Ori
C through the region of the Bar that contains our FOV and outward into the
EON. Not all objects lie in the same direction from $\rm\theta^{1}$ Ori C, so
this two dimensional sketch is accurate only in the sense of the pattern of
irradiation discussed in this section. The smooth line represents the MIF.
O’Dell et al. (2017b) most recently established that inside (closer to
$\rm\theta^{1}$ Ori C) the Bar the MIF is created by EUV from $\rm\theta^{1}$
Ori C. As the MIF tilts toward the observer limb brightening increases the
observed surface brightness until reaching a maximum along the MIF’s
ionization front. Beyond that the MIF begins to flatten, although receiving
enough EUV from $\rm\theta^{1}$ Ori C and later $\rm\theta^{2}$ Ori A to make
the MIF bright. This broad region immediately across the line of maximum tilt
is usually called the Bar (not simply the narrow maximum tilt feature). This
model of the structure of the Bar has been the result of low velocity
resolution spectra that explain the observed ionization changes Mesa-Delgado
et al. (2011), and several high velocity resolution studies employing long-
slit spectra parallel to the Bar made in sequences crossing the Bar (O’Dell et
al., 2017b; O’Dell, 2018). Of course beyond the MIF lies a PDR, that in the
section just outside the visual peak of brightness has produced a cottage-
industry output of insightful infrared and radio studies.

At the outer parts of the shoulder EUV from $\rm\theta^{2}$ Ori A rapidly
begins to be dominant, a condition that increases when moving to greater
distances from $\rm\theta^{1}$ Ori C and $\rm\theta^{2}$ Ori A. The rapid
transition at the shoulder indicates that the shoulder is blocking
$\rm\theta^{1}$ Ori C’s EUV from reaching the MIF further out, with the
consequence that $\rm\theta^{2}$ Ori A becomes the most important star when
looking beyond the vicinity of the Bar’s shoulder (Rubin et al., 2011; O’Dell
et al., 2017b).

In their spectroscopic study of the Huygens Region García-Díaz & Henney (2007)
found an extended low ionization layer of gas northwest of the Bar,
blueshifted with respect to the MIF. Subsequent studies in other emission-
lines Abel et al. (2019) established that this was an ionized layer that in
the vicinity of the Trapezium was about 0.4 pc closer to the observer than
$\rm\theta^{1}$ Ori C (O’Dell et al., 2020). Although called several names in
early studies, the most recent (O’Dell et al., 2020) name is the nearer
ionized layer (NIL), a name indicating that it is closer to the observer than
the MIF that is located on the far side of $\rm\theta^{1}$ Ori C. The source
of the ionization in the center of the Huygens Region is certainly
$\rm\theta^{1}$ Ori C while further out $\rm\theta^{2}$ Ori A probably plays
an important role. One sees the NIL in spectra of the shoulder region. In the
central Huygens Region, the NIL velocity is about 6$\pm$2 km s-1 (O’Dell et
al., 2020) changing to about 0$\pm$2 km s-1 beyond the Bar (O’Dell et al.,
2017b; O’Dell, 2018). The change of velocity in the outer region probably
indicates that the distance of the NIL is influenced by the rising slope of
the MIF.

The average velocity of the NIL emission-lines is similar to the He I 3889 Å
absorption line seen in the Trapezium stars (average 2.1$\pm$0.6 km s-1) and
in $\rm\theta^{2}$ Ori A (two components, at -2.9 km s-1 and 5.2 km s-1). The
2$\rm{}^{3}S$ lower level of the 3889 Å line is a metastable level populated
by recombinations of ionized He and this absorption line most likely arises
from the ionized NIL. Since the 3889 Å absorption line is present in the
spectrum of $\rm\theta^{2}$ Ori A, the NIL must lie on the observer’s side of
the star. Since the NIL is ionization bounded (optically thick to LyC
radiation) (Abel et al., 2019) this means that the NIL permits only
$\rm\theta^{2}$ Ori A’s FUV to pass.

####  5.2.2 Constraints on the location of 203-504

The bright rim of the teardrop shaped 203-504 is clearly pointed towards
$\rm\theta^{2}$ Ori A, as first shown in Bally et al. (2000). Seen in emission
in multiple ions, it shows no evidence of influence by $\rm\theta^{1}$ Ori C
and the EUV of that star must be blocked by intervening gas. These conditions
are satisfied by the proplyd being close to the MIF at the more distant side
of the shoulder, where $\rm\theta^{1}$ Ori C’s radiation is blocked by the
curved PDR that lies underneath. This is shown in Figure 6.

####  5.2.3 Constraints on the location of 203-506

As discussed in Section 4.2.2 203-506 is irradiated only by FUV, although it
is not clear from which nearby star. As shown in Figure 6, this requires it to
be located on the observer’s side of the NIL. This foreground location is
shared by the large number of silhouette-only proplyds seen within the Huygens
Region. These were originally thought to be a product of blocking by the
foreground Veil neutral gas (O’Dell & Wen, 1994), but with the recognition of
the NIL, it now seems the latter is the more likely blocking agent.

The MIF is much higher surface brightness than the NIL, which means that the
fact that 203-506 is seen in silhouette indicates that it is well into the
foreground from the MIF but says nothing about its position relative to the
NIL.

###  5.3 Atomic carbon associated with 203-506.

External photoevaporation is typically identified based on the proplyd
teardrop shaped ionisation front morphology, especially in lines like H
$\alpha$ and Pa $\alpha$ (Bally et al., 2000; Ricci et al., 2008; Kim et al.,
2016; Haworth et al., 2021). To understand the full extent to which external
photoevaporation operates we need to be able to identify it even when there
isn’t a clear teardrop shaped ionisation front. Specifically, when the wind is
weaker, or when a disc is irradiated only by an FUV radiation field. Haworth &
Owen (2020) used PDR-dynamical models and synthetic observations to study the
observational characteristics of winds without ionisation fronts. They found
that CO is typically not expected to be a good tracer of external
photoevaporative winds because it is dissociated near the base of the flow,
making it hard to distinguish from a Keplerian disc. Conversely, Haworth &
Owen (2020) predicted that C i emission could make a good kinematic tracer of
the inner FUV driven parts of winds, however that has not been demonstrated in
practice yet. Haworth et al. (2022) targeted the C i 3P1-3P0 492.16 GHz line
towards proplyds in NGC 1977 with APEX, but obtained no clear detections. The
most likely explanation for this is thought to be low disc masses, but it
could also be due to carbon depletion or beam dilution.

203-506 is the first known instance of an FUV-only externally driven wind from
a protoplanetary disc with no ionisation front (even the weakly irradiated
114-426 has an extended H $\alpha$ ionisation front Miotello et al., 2012). In
the context of the above discussion, the detection of the [Ci] 8727Å line
towards 203-506 illustrated in Figure 2 is the first clear demonstration that
carbon is abundant in the FUV driven part of the wind (see Appendix B for
discussion on the driving mechanism of that line). Although the spectral
resolution of MUSE is inadequate to measure the expected wind velocity ($\sim
4$ km s-1), this makes 203-506 the most promising known target for validating
the kinematic ALMA predictions for C i winds by Haworth & Owen (2020). Also of
note is that more compact [C i] 8727Å emission is detected off-centre of the
YSO and interior to the I-front of 203-504 (see the upper left panel of Figure
4). We estimate [C i] 8727 Å luminosities of $1.1\times 10^{-5}$ L⊙ and
$6.4\times 10^{-7}$ L⊙ for 203-506 and 203-504 respectively.

Goicoechea et al. (priv. communication) have a spatially unresolved [C i]
3P1-3$P_{0}$ detection with ALMA towards 203-506, finding a peak flux density
of 500 mJy/beam at $0.5^{\prime\prime}$ resolution. This makes 203-506 the
best known target for high spatial and spectral resolution follow up with ALMA
to study the kinematics of the FUV heated component of external
photoevaporative winds.

Figure 7: The lines are Ballabio et al. (2023) model predictions of the [O i]
6300Å luminosity due to external photoevaporation as a function of external
FUV radiation field strength and disc radius. The majority of points are data
for ONC proplyds (Ballabio et al., 2023). The square green points are our new
measurements for 203-504 and 203-506.

###  5.4 [O i] 6300Å emission from 203-504 and 203-506

Outside of the extreme conditions in our two targeted proplyds, [O i] 6300 Å
radiation follows collisional excitations from the ground-state $\rm{}^{3}P$
levels to the $\rm{}^{1}D_{2}$ state 2.0 eV above it. This accounts for the
appearance of this line in the MIF and the edge-on view in the Bar. However
Störzer & Hollenbach (1998) showed that [O i] 6300Å can also be emitted
following the FUV driven photodissociation of OH. This was recently
generalised by Ballabio et al. (2023) who demonstrated that theoretical models
of [O i] 6300Å emission by that mechanism are consistent with the [O i]
luminosities from the PDR (i.e. excluding the contribution from the I-front
and jets) of a large number of ONC propylds and predict that the line could be
used to identify external photoevaporation in distant clusters. A continuum
subtracted map of the [O i] 6300Å emission from our FOV is given in the lower
left panel of Figure 4.

One of the issues with studying the [O i] 6300Å line is that it can originate
from multiple sources, including the PDR (as described above), thermally from
internal winds/an ionisation front and as part of high velocity outflows. With
our data we are unable to spectrally distinguish any low and high velocity
components, but do spatially resolve the discs and outflows so can try to
avoid them that way. In the case of 203-504 we estimate the [O i] line
luminosity using an annulus that excludes the bright central component
associated with the outflow and the star itself. For 203-506 we exclude the
bright spot on the northwest side of the disc in the vicinity of the outflow.
Taking this approach we obtain [O i] 6300Å luminosities of $4\pm 2\times
10^{-6}$L⊙ and $4\pm 2\times 10^{-5}$ L⊙ for 203-504 and 203-506 respectively.
We plot these alongside the models and other ONC proplyd data from Ballabio et
al. (2023) in Figure 7. Both discs have [O i] luminosities consistent with the
models. 203-504 has a luminosity similar to the other ONC proplyds. 203-506 is
brighter than many of the proplyds, and at the high end of the observed range
of luminosities found by (Nisini et al., 2018), but agrees well with the model
prediction for a 100 au disc (which is the observed radius of 203-506). That
203-506 is so bright in [O i] 6300Å further supports the interpretation of
Störzer & Hollenbach (1998) and Ballabio et al. (2023) that [O i] 6300Å is
produced in the PDR following the dissociation of OH and not just thermally in
ionised gas or in the shocks of outflows. Zanesse et al. (in prep.) also use
marginally resolved JWST observations to demonstrate that 203-506 has a large
OH abundance.

##  6 Summary and conclusion

We use VLT-MUSE NFM to study two externally photoevaporating discs and their
outflows, situated towards Orion’s SW-Bright Bar. We draw the following main
conclusions from this work.

1) Based on emission line fluxes and ratios, our field is consistent with the
wider irradiation structure of the surrounding ONC. It shows no evidence of
higher ionisation as seen in the vicinity of $\theta^{1}$ Ori C and sits in a
regime where it is possible that both $\theta^{1}$ Ori C and $\theta^{2}$ Ori
A are contributing at a significant level to the radiation field in the
region. Our field contains a $\Delta$ shaped “Delta Feature” which is a
structure in the Bar that exhibits relatively sharp boundaries towards the NE
(in the direction of $\theta^{2}$ Ori A) and towards the NW (in the direction
of $\theta^{1}$ Ori C).

2) The YSO’s, 203-504 and 203-506, have a projected separation of only 660 au,
but are extremely different in the nature of their irradiation. 203-504 is a
classic bright rimmed proplyd with a teardrop shaped ionisation front due to
EUV+FUV irradiation by $\theta^{2}$ Ori A. The fact that 203-506 has no
ionization boundary indicates that it is not irradiated by the EUV from
$\rm\theta^{1}$ Ori C or $\rm\theta^{2}$ Ori A (but is subject to their FUV),
a condition satisfied if 203-506 lies closer to the observer than the near
ionised layer (NIL). The NIL is about 0.4 pc in the foreground of $\theta^{1}$
C (O’Dell et al., 2020) and is of low surface brightness as compared with the
MIF.

3) We detect [Fe ii] 8617Å emission from the outflows HH 519 (emanating from
203-504) and HH 520 (emanating from 203-506). We find that the feature
previously identified as HH 519 is actually the irradiated rim of the Delta
Feature. The new outflow we discover emanating from 203-504, which retains the
name HH 519, is a compact outflow tilted at least 45∘ toward the line of sight
with red- and blue-shifted velocities of around 130 km s-1. The other outflow,
HH 520 emanating from 203-506, was previously considered to be monopolar, but
we demonstrates that it is in fact bipolar, albeit with one component around a
factor 5-10 fainter than the other.

4) The photoevaporative flow around 203-506 is very bright in [C i] 8727 Å
emission. The line is also detected interior to the I-front of 203-504. These
are the first direct spatially resolved detections of atomic carbon emission
in external photoevaporative winds. Atomic carbon is predicted to be a good
kinematic tracer of the FUV driven component of external winds with ALMA.
203-506 is the best known example to test that prediction, but high spectral
resolution is also required.

The pair of discs studied in this paper highlights that the surrounding
ionisation fronts and shielding can affect the radiation field that discs are
exposed to, even when no longer deeply embedded in their natal material, which
has been the focus of recent theoretical work (Concha-Ramírez et al., 2023;
Qiao et al., 2022, 2023; Wilhelm et al., 2023). Further understanding this
interesting phase in disc evolution and environmental impact will require
additional studies of discs at the periphery of H ii regions.

## Acknowledgements

We thank Giacomo Beccari, Fuyan Bian and the ESO support team for their help
in explaining the instrumental artefacts that appear in the H $\alpha$ data.
We thank Bo Reipurth for discussions on the identifier for HH 519, and for
maintaining the Herbig-Haro catalogue.

TJH acknowledges funding from a Royal Society Dorothy Hodgkin Fellowship and
UKRI guaranteed funding for a Horizon Europe ERC consolidator grant
(EP/Y024710/1). JRG thanks the Spanish MCINN for funding support under grant
PID2019-106110GB-I00. CFM and MLA are funded by the European Union (ERC,
WANDA, 101039452). GB is funded by the European Research Council (ERC) under
the European Union’s Horizon 2020 research and innovation programme (Grant
agreement No. 853022, PEVAP). JSK acknowledge support by the National
Aeronautics and Space Administration under Agreement No. 80NSSC21K0593 for the
program "Alien Earths". Views and opinions expressed are however those of the
author(s) only and do not necessarily reflect those of the European Union or
the European Research Council Executive Agency. Neither the European Union nor
the granting authority can be held responsible for them.

## Data Availability

The VLT-MUSE data has a 12 month proprietory period, after which the data will
be publicly available in the ESO archive. We are happy to share in advance of
that date via correspondence.

##  Appendix A Instrumental artefacts

Our reduced dataset included three small dark artefacts in the H  ii region
and near 203-504 in H $\alpha$. They were identified as artefacts (rather than
tiny globules) because in the individual sub-cubes at different inclination
there was a single feature (rather than three) and that was at different
locations. We contacted the ESO helpdesk, who also engaged with the team at
Paranal, to find an explanation, which is summarised below.

In NFM, a dichroic reflects infrared light towards an infrared low order
sensor, which measures the wave front errors invisible to the laser guide star
adaptive optics system. The presence of the dichroic in the light path of MUSE
when used in the NFM-AO, has some effect on the quality of the signal. The
most noticeable effect of the dichroic is the introduction of a few small
regions (spots) of reduced sensitivity, varying from a few up to 10 per cent
compared to the median value. This is likely triggered by the presence of
impurities or small particles of dust on the dichroic itself. Using flat
fields and images at different rotations efficiently removes such effects from
the combined cubes but sporadically, they can be visible. This is the most
likely suggested cause of the instrumental artefacts in our dataset.

##  Appendix B Explaining the [C i] 8727 Å Emission

[C i] emission is seen in spectra of the Huygens region in the form of the
9824.31 Å, 9850.34 Å and 8727.126 Å lines (Esteban et al., 2004). The longer
wavelength lines are forbidden transitions from the 2$\rm s^{2}$ 2$\rm p^{2}$
$\rm{}^{1}D_{2}$ level at 1.26 eV above the ground state. The shorter
wavelength line is a forbidden transition from the 2$\rm s^{2}$ 2$\rm p^{2}$
$\rm{}^{1}S_{0}$ level 2.7 eV above the ground state to the upper level of the
longer wavelength lines.

With an ionisation potential of 11.26 eV, atomic carbon is ionised by FUV in
much of the PDR. Escalante et al. (1991) argue that recombinations populate
the state producing the 8727 Å line, as well as the state giving rise to the
9850 Å and 9824 Å lines. In addition, collisions are also a possible origin of
the [C i] 8727 Å emission, however, collisional excitation in the $\sim 9000$
K gas of the Huygens region should lead to 9850 Å and 9824 Å lines being much
stronger than the [C i] 8727 Å emission, which Esteban et al. (2004) show is
not the case. Specifically, Esteban et al. (2004) find the 9824 Å / 8727 Å
flux ratio is of order unity and the 9850 Å / 8727 Å flux ratio is $\sim 5$,
which is consistent with recombination excitation in a PDR (Escalante et al.,
1991). We hence conclude that the most likely source of the 8727 Å emission is
due to photoionizations by FUV followed by recombinations rather than
collisional excitation.

##  Appendix C Explaining the [O i] 8446 Å Emission

The surprising presence of [O i] 8446 Å lines in the spectra of gaseous
nebulae was first explained as the result of fluorescence by Ly$\beta$ by
Bowen (1928), with the details of the process being explained more recently in
Section 4.7 and Figure 4.7 of Osterbrock & Ferland (2006). Ly$\beta$ (1025.72
Å) has a near coincidence in wavelength with a ground state line of O i at
1025.76 Å, driving excitation of the upper state of the neutral oxygen,
followed by spontaneous decay. Production of the line requires both the
presence of neutral oxygen and Ly$\beta$ photons from ionized hydrogen. The
Ly$\beta$ photons could come from co-existing ionized hydrogen interacting
with neutral oxygen produced by the mechanism of charge exchange, where
ionized oxygen gain an electron from the highly ionized hydrogen gas. An
alternative source of Ly$\beta$ is emission from the ionized hydrogen zone
that begins at the ionization boundary.

The ionization balance in the Huygens Region is very dependent on the distance
from $\rm\theta^{1}$ Ori C, as shown in the multi-sample study of O’Dell &
Harris (2010). Because of problems of correcting for night-sky [O i] 6300 Å
the pattern of appearance of [O i] could not be traced, but the pattern of
increasing relative strength of low ionization lines is shown in their Figure
11 ([S ii]) and Figure 12 ([N ii]). This pattern of increasing relative
strength of [S ii] and [N ii] leads one to expect an increasing relative
strength of [O i], which explains the larger relative strengths of [O i] 8446
Å in the outer regions in Table 4.

Figure 8: The ratio of [O i] 8446 Å to [O i] 6300 Å.The strength of [O i] 8446
Å is driven by resonance with Lyman $\beta$ , c. f. Osterbrock & Ferland
(2006) and that of [O i] 6300 Å by collisional excitation (except within the
external photoevaporative winds) as explained in section 5.4. The variations
are largely caused by velocity differences of the emitter and absorber ions.

The relative velocities of the Ly$\beta$ emitting (hydrogen) and Ly$\beta$
absorbing (oxygen) gas plays a role in [O i] 8446 Å emission.If both oxygen
and hydrogen have the same average velocity, the peak of the oxygen line
profile is 11.7 km s-1 to the red from the peak of the Ly$\beta$ line.At 9000
K, the full width at half maximum intensity is 20.4 km s-1 for hydrogen and
5.1 km s-1 for oxygen. The imperfect alignment of the two lines (the oxygen
line peak is about at the half-power point of the hydrogen profile) means that
the efficiency of the production of [O i] is sensitive to their relative
velocities. For example, shifting the oxygen to a lower velocity difference
would increase the strength of the 8446 Å line. Table 5 of O’Dell et al.
(2020) summarizes the velocities of difference ions in the MIF in the central
Huygens Region, giving 27$\pm$2 km s-1for [O i] and 17$\pm$2 km s-1 for
hydrogen Balmer lines. This means that if the Ly$\beta$ photons pumping [O i]
8446 Å arise from a different zone, the separation of the peaks of the two
lines has increased by 10 km s-1 and the enhancement of [O i] 8446 Å would be
reduced, as compared with a condition where the Ly$\beta$ source and the
neutral oxygen have the same velocity.

The velocity difference in the central Huygens Region is due to
photoevaporative flow away from the hydrogen ionization front that is being
viewed face-on. If the flow velocity is less in this outer region where the
EUV is weaker, then the [O i] 8446 Å emissivity would be increased, accounting
for the enhancement seen in Table 4.

Differences of relative velocities are likely to play a role within the Delta
Feature. In Figure 8 we show the ratio of the [O i] 8446Å to the [O i] 6300Å
line. Ignoring the photoevaporating discs (where processes such as OH
dissociation determine the 6300 Å line flux) in the surrounding medium we see
the results of differences of velocity of the emitting and absorbing gas, with
the more intense [O i] 8446 Å arising from greater overlap of the Ly$\beta$
and neutral oxygen line profiles while the [O i] 6300 emission arises from
collisional excitation by the electron gas, independent of the gas velocity.
At the well defined NE and NW boundaries of the Delta Feature abrupt changes
of relative velocity are to be expected, as shown in high velocity resolution
studies of the Bar (O’Dell et al., 2017a), if this feature is a discrete
physical structure. The slowly changing values of the ratio along the south
boundary of the Delta Feature are consistent with the idea that the boundary
there is more gradual.

## References

  * ALMA Partnership et al. (2015) ALMA Partnership et al., 2015, ApJ, 808, L3
  * Abel et al. (2019) Abel N. P., Ferland G. J., O’Dell C. R., 2019, ApJ, 881, 130
  * Abt et al. (1991) Abt H. A., Wang R., Cardona O., 1991, ApJ, 367, 155
  * Adams et al. (2004) Adams F. C., Hollenbach D., Laughlin G., Gorti U., 2004, ApJ, 611, 360
  * Agra-Amboage et al. (2011) Agra-Amboage V., Dougados C., Cabrit S., Reunanen J., 2011, A&A, 532, A59
  * Andree-Labsch et al. (2017) Andree-Labsch S., Ossenkopf-Okada V., Röllig M., 2017, A&A, 598, A2
  * Andrews et al. (2018) Andrews S. M., et al., 2018, ApJ, 869, L41
  * Bacon et al. (2010) Bacon R., et al., 2010, in McLean I. S., Ramsay S. K., Takami H., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 7735, Ground-based and Airborne Instrumentation for Astronomy III. p. 773508 (arXiv:2211.16795), doi:10.1117/12.856027
  * Baldwin et al. (1991) Baldwin J. A., Ferland G. J., Martin P. G., Corbin M. R., Cota S. A., Peterson B. M., Slettebak A., 1991, ApJ, 374, 580
  * Ballabio et al. (2023) Ballabio G., Haworth T. J., Henney W. J., 2023, MNRAS, 518, 5563
  * Bally et al. (2000) Bally J., O’Dell C. R., McCaughrean M. J., 2000, AJ, 119, 2919
  * Berné et al. (2022) Berné O., et al., 2022, PASP, 134, 054301
  * Blagrave et al. (2007) Blagrave K. P. M., Martin P. G., Rubin R. H., Dufour R. J., Baldwin J. A., Hester J. J., Walter D. K., 2007, ApJ, 655, 299
  * Bowen (1928) Bowen I. S., 1928, ApJ, 67, 1
  * Boyden & Eisner (2020) Boyden R. D., Eisner J. A., 2020, ApJ, 894, 74
  * Boyden & Eisner (2023) Boyden R. D., Eisner J. A., 2023, ApJ, 947, 7
  * Bregman et al. (1989) Bregman J. D., Allamandola L. J., Tielens A. G. G. M., Geballe T. R., Witteborn F. C., 1989, ApJ, 344, 791
  * Cardelli et al. (1989) Cardelli J. A., Clayton G. C., Mathis J. S., 1989, ApJ, 345, 245
  * Champion et al. (2017) Champion J., Berné O., Vicente S., Kamp I., Le Petit F., Gusdorf A., Joblin C., Goicoechea J. R., 2017, A&A, 604, A69
  * Chandar (2009) Chandar R., 2009, Ap&SS, 324, 315
  * Clarke & Owen (2015) Clarke C. J., Owen J. E., 2015, MNRAS, 446, 2944
  * Concha-Ramírez et al. (2023) Concha-Ramírez F., Wilhelm M. J. C., Portegies Zwart S., 2023, MNRAS, 520, 6159
  * Cuadrado et al. (2015) Cuadrado S., Goicoechea J. R., Pilleri P., Cernicharo J., Fuente A., Joblin C., 2015, A&A, 575, A82
  * Cuello et al. (2019) Cuello N., et al., 2019, MNRAS, 483, 4114
  * Cuello et al. (2020) Cuello N., et al., 2020, MNRAS, 491, 504
  * Cuello et al. (2023) Cuello N., Ménard F., Price D. J., 2023, European Physical Journal Plus, 138, 11
  * Dale & Bonnell (2011) Dale J. E., Bonnell I., 2011, MNRAS, 414, 321
  * Dale et al. (2012) Dale J. E., Ercolano B., Bonnell I. A., 2012, MNRAS, 424, 377
  * Dale et al. (2014) Dale J. E., Ngoumou J., Ercolano B., Bonnell I. A., 2014, MNRAS, 442, 694
  * Dipierro et al. (2015) Dipierro G., Price D., Laibe G., Hirsh K., Cerioli A., Lodato G., 2015, MNRAS, 453, L73
  * Eisner et al. (2018) Eisner J. A., et al., 2018, ApJ, 860, 77
  * Escalante et al. (1991) Escalante V., Sternberg A., Dalgarno A., 1991, ApJ, 375, 630
  * Esteban et al. (2004) Esteban C., Peimbert M., García-Rojas J., Ruiz M. T., Peimbert A., Rodríguez M., 2004, MNRAS, 355, 229
  * Facchini et al. (2016) Facchini S., Clarke C. J., Bisbas T. G., 2016, MNRAS, 457, 3593
  * Fall et al. (2010) Fall S. M., Krumholz M. R., Matzner C. D., 2010, ApJ, 710, L142
  * Fang et al. (2016) Fang M., Kim J. S., Pascucci I., Apai D., Manara C. F., 2016, ApJ, 833, L16
  * Fang et al. (2021) Fang M., Kim J. S., Pascucci I., Apai D., 2021, ApJ, 908, 49
  * Ferland et al. (1998) Ferland G. J., Korista K. T., Verner D. A., Ferguson J. W., Kingdon J. B., Verner E. M., 1998, PASP, 110, 761
  * Ferland et al. (2013) Ferland G. J., et al., 2013, Rev. Mex. Astron. Astrofis., 49, 137
  * Ferland et al. (2017) Ferland G. J., et al., 2017, Rev. Mex. Astron. Astrofis., 53, 385
  * Freudling et al. (2013) Freudling W., Romaniello M., Bramich D. M., Ballester P., Forchi V., García-Dabló C. E., Moehler S., Neeser M. J., 2013, A&A, 559, A96
  * Fukui et al. (2021) Fukui Y., Habe A., Inoue T., Enokiya R., Tachihara K., 2021, PASJ, 73, S1
  * Gaia Collaboration et al. (2016) Gaia Collaboration et al., 2016, A&A, 595, A1
  * Gaia Collaboration et al. (2018) Gaia Collaboration et al., 2018, A&A, 616, A1
  * Gaia Collaboration et al. (2022) Gaia Collaboration et al., 2022, arXiv e-prints, p. arXiv:2208.00211
  * García-Díaz & Henney (2007) García-Díaz M. T., Henney W. J., 2007, AJ, 133, 952
  * Garufi et al. (2022) Garufi A., et al., 2022, A&A, 658, A104
  * Gaudi et al. (2021) Gaudi B. S., Meyer M., Christiansen J., 2021, in Madhusudhan N., ed., , ExoFrontiers; Big Questions in Exoplanetary Science. pp 2–1, doi:10.1088/2514-3433/abfa8fch2
  * Giannini et al. (2013) Giannini T., et al., 2013, ApJ, 778, 71
  * Ginski et al. (2021) Ginski C., et al., 2021, ApJ, 908, L25
  * Goicoechea et al. (2016) Goicoechea J. R., et al., 2016, Nature, 537, 207
  * Greve et al. (1994) Greve A., Castles J., McKeith C. D., 1994, A&A, 284, 919
  * Gritschneder et al. (2009) Gritschneder M., Naab T., Walch S., Burkert A., Heitsch F., 2009, ApJ, 694, L26
  * Gritschneder et al. (2010) Gritschneder M., Burkert A., Naab T., Walch S., 2010, ApJ, 723, 971
  * Großschedl et al. (2018) Großschedl J. E., et al., 2018, A&A, 619, A106
  * Grudić et al. (2021) Grudić M. Y., Guszejnov D., Hopkins P. F., Offner S. S. R., Faucher-Giguère C.-A., 2021, MNRAS, 506, 2199
  * Grudić et al. (2022) Grudić M. Y., Guszejnov D., Offner S. S. R., Rosen A. L., Raju A. N., Faucher-Giguère C.-A., Hopkins P. F., 2022, MNRAS, 512, 216
  * Gupta et al. (2023) Gupta A., et al., 2023, A&A, 670, L8
  * Guszejnov et al. (2022) Guszejnov D., Grudić M. Y., Offner S. S. R., Faucher-Giguère C.-A., Hopkins P. F., Rosen A. L., 2022, MNRAS, 515, 4929
  * Habart et al. (2023) Habart E., et al., 2023, A&A, 673, A149
  * Haffert et al. (2019) Haffert S. Y., Bohn A. J., de Boer J., Snellen I. A. G., Brinchmann J., Girard J. H., Keller C. U., Bacon R., 2019, Nature Astronomy, 3, 749
  * Hartigan & Morse (2007) Hartigan P., Morse J., 2007, ApJ, 660, 426
  * Hartigan et al. (2001) Hartigan P., Morse J. A., Reipurth B., Heathcote S., Bally J., 2001, ApJ, 559, L157
  * Hartigan et al. (2004) Hartigan P., Raymond J., Pierson R., 2004, ApJ, 614, L69
  * Hartigan et al. (2005) Hartigan P., Heathcote S., Morse J. A., Reipurth B., Bally J., 2005, AJ, 130, 2197
  * Haworth & Clarke (2019) Haworth T. J., Clarke C. J., 2019, MNRAS, 485, 3895
  * Haworth & Owen (2020) Haworth T. J., Owen J. E., 2020, MNRAS, 492, 5030
  * Haworth et al. (2018) Haworth T. J., Clarke C. J., Rahman W., Winter A. J., Facchini S., 2018, MNRAS, 481, 452
  * Haworth et al. (2021) Haworth T. J., Kim J. S., Winter A. J., Hines D. C., Clarke C. J., Sellek A. D., Ballabio G., Stapelfeldt K. R., 2021, MNRAS, 501, 3502
  * Haworth et al. (2022) Haworth T. J., et al., 2022, MNRAS, 512, 2594
  * Henney & O’Dell (1999) Henney W. J., O’Dell C. R., 1999, AJ, 118, 2350
  * Henney et al. (2005) Henney W. J., Arthur S. J., García-Díaz M. T., 2005, ApJ, 627, 813
  * Henroteau & Henderson (1920) Henroteau F., Henderson J. P., 1920, Publications of the Dominion Observatory Ottawa, 5, 1
  * Huang et al. (2022) Huang J., et al., 2022, ApJ, 930, 171
  * Jansen et al. (1995) Jansen D. J., Spaans M., Hogerheijde M. R., van Dishoeck E. F., 1995, A&A, 303, 541
  * Jeffries (2007) Jeffries R. D., 2007, MNRAS, 376, 1109
  * Johnstone et al. (1998) Johnstone D., Hollenbach D., Bally J., 1998, ApJ, 499, 758
  * Keppler et al. (2018) Keppler M., et al., 2018, A&A, 617, A44
  * Kim et al. (2016) Kim J. S., Clarke C. J., Fang M., Facchini S., 2016, ApJ, 826, L15
  * Kirwan et al. (2023) Kirwan A., et al., 2023, A&A, 673, A166
  * Kounkel et al. (2017) Kounkel M., et al., 2017, ApJ, 834, 142
  * Kounkel et al. (2018) Kounkel M., et al., 2018, AJ, 156, 84
  * Kratter & Lodato (2016) Kratter K., Lodato G., 2016, ARA&A, 54, 271
  * Kurucz (1979) Kurucz R. L., 1979, ApJS, 40, 1
  * Lada & Lada (2003) Lada C. J., Lada E. A., 2003, ARA&A, 41, 57
  * Luridiana et al. (2015) Luridiana V., Morisset C., Shaw R. A., 2015, A&A, 573, A42
  * Mann et al. (2014) Mann R. K., et al., 2014, ApJ, 784, 82
  * McLeod et al. (2016) McLeod A. F., Weilbacher P. M., Ginsburg A., Dale J. E., Ramsay S., Testi L., 2016, MNRAS, 455, 4057
  * Mesa-Delgado et al. (2011) Mesa-Delgado A., Núñez-Díaz M., Esteban C., López-Martín L., García-Rojas J., 2011, MNRAS, 417, 420
  * Miotello et al. (2012) Miotello A., Robberto M., Potenza M. A. C., Ricci L., 2012, ApJ, 757, 78
  * Miotello et al. (2021) Miotello A., Rosotti G., Ansdell M., Facchini S., Manara C. F., Williams J. P., Bruderer S., 2021, A&A, 651, A48
  * Nisini et al. (2018) Nisini B., Antoniucci S., Alcalá J. M., Giannini T., Manara C. F., Natta A., Fedele D., Biazzo K., 2018, A&A, 609, A87
  * O’Dell (2018) O’Dell C. R., 2018, MNRAS, 478, 1017
  * O’Dell & Harris (2010) O’Dell C. R., Harris J. A., 2010, AJ, 140, 985
  * O’Dell & Wen (1994) O’Dell C. R., Wen Z., 1994, ApJ, 436, 194
  * O’Dell et al. (1993) O’Dell C. R., Wen Z., Hu X., 1993, ApJ, 410, 696
  * O’Dell et al. (2009) O’Dell C. R., Henney W. J., Abel N. P., Ferland G. J., Arthur S. J., 2009, AJ, 137, 367
  * O’Dell et al. (2017a) O’Dell C. R., Ferland G. J., Peimbert M., 2017a, MNRAS, 464, 4835
  * O’Dell et al. (2017b) O’Dell C. R., Kollatschny W., Ferland G. J., 2017b, ApJ, 837, 151
  * O’Dell et al. (2020) O’Dell C. R., Abel N. P., Ferland G. J., 2020, ApJ, 891, 46
  * Öberg et al. (2021) Öberg K. I., et al., 2021, ApJS, 257, 1
  * O’dell & Wong (1996) O’dell C. R., Wong K., 1996, AJ, 111, 846
  * Osterbrock & Ferland (2006) Osterbrock D. E., Ferland G. J., 2006, Astrophysics of gaseous nebulae and active galactic nuclei
  * Osterbrock et al. (1992) Osterbrock D. E., Tran H. D., Veilleux S., 1992, ApJ, 389, 305
  * Pesenti et al. (2003) Pesenti N., Dougados C., Cabrit S., O’Brien D., Garcia P., Ferreira J., 2003, A&A, 410, 155
  * Pfalzner (2013) Pfalzner S., 2013, A&A, 549, A82
  * Pfalzner et al. (2005) Pfalzner S., Umbreit S., Henning T., 2005, ApJ, 629, 526
  * Pfalzner et al. (2018) Pfalzner S., Bhandare A., Vincke K., Lacerda P., 2018, ApJ, 863, 45
  * Pineda et al. (2020) Pineda J. E., Segura-Cox D., Caselli P., Cunningham N., Zhao B., Schmiedeke A., Maureira M. J., Neri R., 2020, Nature Astronomy, 4, 1158
  * Pinte et al. (2018) Pinte C., et al., 2018, ApJ, 860, L13
  * Pinte et al. (2020) Pinte C., et al., 2020, ApJ, 890, L9
  * Pinte et al. (2022) Pinte C., Teague R., Flaherty K., Hall C., Facchini S., Casassus S., 2022, arXiv e-prints, p. arXiv:2203.09528
  * Portegies Zwart et al. (2010) Portegies Zwart S. F., McMillan S. L. W., Gieles M., 2010, ARA&A, 48, 431
  * Qiao et al. (2022) Qiao L., Haworth T. J., Sellek A. D., Ali A. A., 2022, MNRAS, 512, 3788
  * Qiao et al. (2023) Qiao L., Coleman G. A. L., Haworth T. J., 2023, MNRAS, 522, 1939
  * Reiter & Parker (2022) Reiter M., Parker R. J., 2022, European Physical Journal Plus, 137, 1071
  * Reiter & Smith (2013) Reiter M., Smith N., 2013, MNRAS, 433, 2226
  * Reiter et al. (2015) Reiter M., Smith N., Kiminki M. M., Bally J., 2015, MNRAS, 450, 564
  * Reiter et al. (2016) Reiter M., Smith N., Bally J., 2016, MNRAS, 463, 4344
  * Reiter et al. (2017) Reiter M., Kiminki M. M., Smith N., Bally J., 2017, MNRAS, 470, 4671
  * Reiter et al. (2019) Reiter M., McLeod A. F., Klaassen P. D., Guzmán A. E., Dale J. E., Mottram J. C., Garay G., 2019, MNRAS, 490, 2056
  * Reiter et al. (2022) Reiter M., Morse J. A., Smith N., Haworth T. J., Kuhn M. A., Klaassen P. D., 2022, MNRAS, 517, 5382
  * Ricci et al. (2008) Ricci L., Robberto M., Soderblom D. R., 2008, AJ, 136, 2136
  * Richling & Yorke (2000) Richling S., Yorke H. W., 2000, ApJ, 539, 258
  * Rubin et al. (2011) Rubin R. H., Simpson J. P., O’Dell C. R., McNabb I. A., Colgan S. W. J., Zhuge S. Y., Ferland G. J., Hidalgo S. A., 2011, MNRAS, 410, 1320
  * Schulz et al. (2006) Schulz N. S., Testa P., Huenemoerder D. P., Ishibashi K., Canizares C. R., 2006, ApJ, 653, 636
  * Sicilia-Aguilar et al. (2005) Sicilia-Aguilar A., et al., 2005, AJ, 129, 363
  * Storey & Hummer (1995) Storey P. J., Hummer D. G., 1995, MNRAS, 272, 41
  * Störzer & Hollenbach (1998) Störzer H., Hollenbach D., 1998, ApJ, 502, L71
  * Tang et al. (2012) Tang Y. W., Guilloteau S., Piétu V., Dutrey A., Ohashi N., Ho P. T. P., 2012, A&A, 547, A84
  * Tauber et al. (1994) Tauber J. A., Tielens A. G. G. M., Meixner M., Goldsmith P. F., 1994, ApJ, 422, 136
  * Teague et al. (2018) Teague R., Bae J., Bergin E. A., Birnstiel T., Foreman-Mackey D., 2018, ApJ, 860, L12
  * Valdivia-Mena et al. (2022) Valdivia-Mena M. T., et al., 2022, A&A, 667, A12
  * Walch et al. (2012) Walch S. K., Whitworth A. P., Bisbas T., Wünsch R., Hubber D., 2012, MNRAS, 427, 625
  * Walmsley et al. (2000) Walmsley C. M., Natta A., Oliva E., Testi L., 2000, A&A, 364, 301
  * Weilbacher et al. (2015) Weilbacher P. M., et al., 2015, A&A, 582, A114
  * Wen & O’Dell (1995) Wen Z., O’Dell C. R., 1995, ApJ, 438, 784
  * Wilhelm et al. (2023) Wilhelm M. J. C., Portegies Zwart S., Cournoyer-Cloutier C., Lewis S. C., Polak B., Tran A., Mac Low M.-M., 2023, MNRAS, 520, 5331
  * Winter & Haworth (2022) Winter A. J., Haworth T. J., 2022, European Physical Journal Plus, 137, 1132
  * Winter et al. (2019) Winter A. J., Clarke C. J., Rosotti G. P., Hacar A., Alexander R., 2019, MNRAS, 490, 5478
  * Yen et al. (2019) Yen H.-W., Gu P.-G., Hirano N., Koch P. M., Lee C.-F., Liu H. B., Takakuwa S., 2019, ApJ, 880, 69
  * Zeidler (2019) Zeidler P., 2019, MUSEpack, Zenodo, doi:10.5281/zenodo.3433996
  * Zeidler et al. (2019) Zeidler P., Nota A., Sabbi E., Luljak P., McLeod A. F., Grebel E. K., Pasquali A., Tosi M., 2019, AJ, 158, 201
  * Zhang et al. (2018) Zhang S., et al., 2018, ApJ, 869, L47
  * Zhang et al. (2023) Zhang S., Kalscheur M., Long F., Zhang K., Long D. E., Bergin E. A., Zhu Z., Trapman L., 2023, arXiv e-prints, p. arXiv:2305.03862
  * Zucker et al. (2022) Zucker C., et al., 2022, Nature, 601, 334