# The Young Stellar Population, Distance, and Cloud–Cloud Collision Induced
Star Formation Scenario of the Trifid Nebula

Venu M. Kalari<EMAIL_ADDRESS>Gemini Observatory/NSF’s NOIRLab,
Casilla 603, La Serena, Chile  Departamento de Astronomia, Universidad de
Chile, Casilla 36-D, Santiago, Chile

###### Abstract

The Trifid Nebula is a young, nearby star-forming region where star formation
is proposed to have been triggered by cloud-cloud collision (CCC), based on
observations of molecular clouds. It offers a unique opportunity to test
whether the CCC hypothesis is supported by the spatial distribution and star
formation chronology of young stars.

We present the first study of the optically visible pre-main sequence (PMS)
population of the region using riH$\alpha$ imaging and Gaia astrometry.
Combined with an analysis of young stellar objects (YSOs) using infrared
imaging, we capture the spatial distribution and star formation chronology of
the young stellar population. From the analysis, 15 Flat/Class I YSOs, 46
Class II YSOs, and 41 accreting PMS stars are identified (diskless/non-
accreting sources are not included in the analysis). The distance based on
Gaia parallaxes is $\sim$1250 pc, significantly closer than previously
reported. The Class II YSOs and PMS stars ($\sim$1.5 Myr old) are spread
toward the edge of the molecular clouds. They are slightly younger than the
estimated crossing time of $\sim$2.7 Myr and closer to the estimated dynamical
age $\sim$0.85 Myr. Younger Class I YSOs are more concentrated spatially.
There exists a cavity devoid of young stars where the two clouds overlap. This
evidence suggests that the current generation of stars formed after the
collision of two clouds $\sim$1 Myr ago, and this result can be corroborated
using future spectroscopic studies.

H II regions: Trifid Nebula — stars: formation — stars: pre-main sequence

††facilities: Spitzer,Herschel, VLT(VST), Gaia††software: APLpy (Robitaille &
Bressert, 2012), FILFinder (Koch & Rosolowsky, 2015)

##  1 Introduction

Figure 1: The Trifid nebula as an rgb image from VPHAS+ $i$, H$\alpha$, and
$r$ images, respectively. The area shown is centred on right ascension (J2000)
18h02m26s8 and declination (J2000) $-$23$\arcdeg$00$\arcmin$12$\arcsec$ with a
field of view 12$\arcmin$2 $\times$ 14$\arcmin$9, and covers the area studied
in this work. North is up and east to the left. The central ionising O7.5V
star in the HD 164492 system is visible as the bright central star.
Nebulosity, varying on small spatial scales is also seen.

There remain many open questions in the field of star formation. In the
current paradigm, star formation is a dynamical process, occurring primarily
due to gravitational instabilities seeded by turbulence in molecular clouds
(McKee & Ostriker, 2007). A spectacular trigger for gravitational
instabilities and the eventual collapse of molecular clouds into stars is
cloud-cloud collision (CCC). It is has been argued that the collision of one
or more molecular clouds can potentially cause over-densities, leading to
gravitational collapse and subsequent star formation (Loren, 1976). Recent
observations at sub-mm wavelengths have presented evidence that CCC is a
viable mechanism to form stars and young clusters across the Milky Way and in
external galaxies (Fukui et al., 2020).

The young star-forming region Trifid Nebula (or M 20, NGC 6514) is a proposed
site of CCC induced star formation (Torii et al., 2011, 2017), based on the
systemic velocities and masses of two molecular clouds in the region. These
two clouds are thought to have collided $\sim$1 Myr ago, producing the young
stellar population observed in the region. M 20 is one of the youngest sites
($\sim$0.3–1 Myr; Rho et al. 2008) proposed to have formed via CCC. It is
thought to lie between 1.6–2.7 kpc away (Rho et al., 2008; Cambrésy et al.,
2011), either in the Sagittarius or Scutum arm. This combination of proximity
and age allows one to observe and characterise the young low-mass stellar
population of the region, a possibility not afforded in studies of more
distant sites of CCC (Fukui et al., 2020). It is therefore an ideal avenue to
gauge whether the young stellar population is supportive of an induced star-
formation scenario.

The stellar population of the region has not been studied in the optical. An
O7.5V star, HD 164492A (with a spectroscopic age of 0.6 Myr; Petit et al.
2019) ionises the optically observed nebula (see Fig. 1), and lies south west
of the crux of the three bisecting dust lanes. The extinction law in the
region is anomalous at $R_{V}=5.5$ (Cambrésy et al., 2011). The combination of
nebulousity, extinction, and a line-of-sight towards the Galactic plane have
hindered optical observations. Although there have been no optical studies of
the young stellar population in the region, multiple studies at radio, sub-mm,
mid-infrared (mIR), and near-infrared (nIR) have been conducted (Yusef-Zadeh
et al., 2000; Lefloch et al., 2002; Rho et al., 2001, 2006, 2008; Feigelson et
al., 2013). These studies have revealed a rich young stellar population at a
very nascent stage of star formation ($\lesssim$1 Myr old); a population of
even younger star-forming cores (Tapia et al., 2018), along with evidence for
photo-ionisation affecting on-going star formation (Yusef-Zadeh et al., 2005)

The different wavelengths trace the varying stages of on-going star formation
in the region; from the molecular clouds visible at sub-mm and radio
wavelengths, to the young protostars in the infrared (IR). Studies of the
optically visible pre-main sequence (PMS) stars and the stellar cluster would
complete the circle, helping to piece together the star formation picture.
With a combination of deep optical photometry, and astrometry from the Gaia
space telescope, it is now possible to both observe the low-mass optically
visible PMS stars and disentangle them from line-of-sight of contaminants.
Thus, when combined with archival IR data, we are now armed with observations
to characterise the young stellar population of the Trifid Nebula, and examine
how they relate to the proposed CCC triggered star formation scenario.

In this paper we identify and characterise the young stellar population in the
Trifid Nebula using IR and optical observations, combined with Gaia astrometry
to investigate whether the triggered star formation scenario proposed in the
literature can explain the observed spatial distribution and star formation
chronology of the young stellar population. This paper is organised thus; in
Section 2, a description of the data used in this work is presented. Sections
3 and 4 describe the identification and classification of the IR young stellar
population, and optically visible PMS stars respectively. In Section 5 a
discussion on the star-formation scenario in the Trifid nebula is presented.

##  2 Data

Figure 2: Observed ($r-i$) vs. $r$ colour magnitude diagram of all sources
studied in this work. A reddening vector of $A_{V}$=1 mag, $R_{V}$=5.5, is
shown in the top right. Bressan et al. (2012) isochrones of 1 (solid) and 10
Myr (dashed) at a distance of 1250 pc and $A_{V}$ = 1.3 mag are shown.

###  2.1 Area of study

To identify and characterise the young stellar population in the star-forming
region M 20, we define our area of study as all sources falling within a
bounding box of 12$\arcmin$2 $\times$ 14$\arcmin$9 region centred on right
ascension 18h02m26s8 and declination $-$23$\arcdeg$00$\arcmin$12$\arcsec$ (see
Fig. 1). The area definition is slightly smaller than previous studies at nIR
or X-ray wavelengths (Rho et al., 2001, 2006; Feigelson et al., 2013), and
completely includes only the central Trifid nebula, and the blister region
located immediately above it. Because M 20 lies directly in our line of sight
towards the Galactic plane, multiple populations in further arms will be
confused with optical members. The optical nebulosity as seen in M 20 can
prevent the detection of more distant contaminants in the plane, providing an
upper blue envelope in the colour-magnitude diagram. In addition, the area
chosen adequately covers the region where CCC and subsequent star formation
are thought to have occurred (Torii et al., 2011). Choosing a slightly larger
area of study while potentially enabling a more detailed understanding of the
extent of the region, will include more distant non-members which cannot be
separated easily, justifying our choice of boundaries.

###  2.2 VPHAS+ ugriH$\alpha$ imaging

The VPHAS+ (VST/Omegacam Photometric H$\alpha$ Survey of the Southern galactic
plane and bulge) survey observed the Trifid Nebula in ugriH$\alpha$ filters,
and is comprehensively outlined in the survey paper by Drew et al. (2014). The
data used in this publication is part of the VPHAS+ third data release. The
reader is referred to the release
document111http://www.eso.org/rm/api/v1/public/releaseDescriptions/106 for
details of the data reduction procedure. A rgb image created using the final
reduced images is shown in Fig.1 to illustrate the quality of the imaging and
the area studied in this work. Clearly visible are the three central dust
lanes, and the photoionized regions towards the south.

Aperture photometry was performed on these reduced images using a dedicated
pipeline run by the Cambridge Astronomy Survey Unit. Attention is drawn to the
dedicated software nebuliser used in the photometric pipeline to track
variable background emission on a sliding scale around 15$\arcsec$. Objects
that suffer from improper background subtraction due to rapidly varying
nebulousity on smaller scales are flagged using a curve-of-growth analysis as
possibly non-stellar. The final catalogue available from the archive is
calibrated in the AB magnitude system. A global Vega magnitude calibration is
achieved across all filters by cross matching sources overlapping with the
VPHAS+ DR2 data , which was globally calibrated on the Vega magnitude scale.
The global shifts are 1.084, $-$0.055, 0.138, 0.461, and 0.316 in
ugriH$\alpha$ filters respectively. For the remainder of the study, the
magnitudes reported are in the Vega magnitude scale.

From the publicly available catalogue, we select all sources falling within
the area of interest. We apply further selection criteria of (i) $r>$13 to
remove saturated sources, (ii) photometric uncertainties in $ri$ less than 0.1
mag, and in H$\alpha$ less than 0.15 mag to keep low propagated errors, (iii)
sources classified as stellar, or star-like in broadband $ri$ photometry to
remove extended objects. For H$\alpha$ photometry, mild degradation of the
point-spread function may lead to spurious classifications as extended
objects. To avoid this and not reject acceptable H$\alpha$ measurements, the
classification in H$\alpha$ is relaxed to extended sources. Instead, a further
visual check using unsharp masking on selected sources is imposed to eliminate
any sources affected by sharply varying nebulousity (see Kalari et al. 2015),
(iv) sources with average photometric confidence greater than 90 to remove
sources falling in the chip gaps. Note that we do not impose any quality
criteria on $ug$ photometry. The cleaned sample contains 9 227 unique sources
with $ri$H$\alpha$ photometry meeting the quality criteria. The colour
magnitude diagram generated using the data is shown in Fig. 2. The saturation
limit of VPHAS+ photometry means that the bright, early type members of the
Trifid nebula are missing from our analysis. Unlike other nearby star-forming
regions near the Galactic plane, there is no upper blue envelope due to the
nebulosity. More distant Galactic contaminants are visible, particularly at
the edges of the nebula.

Finally, a result of our selection criteria is that stars near the nebulous
centre, and along the dust lanes are preferentially removed, as they lie along
regions of very rapidly spatially varying nebulousity. Optical photometry from
the Gaia telescope also shows similar spatial density to the VPHAS+
observations, suggesting that observations at the optical wavelength towards
the central region suffer greatly from the spatial varying nebulousity. The
consequences of this spatial density variation in the optical on our results
is discussed in Section 4.6.

###  2.3 Gaia EDR 3 astrometry

Astrometric data from the Gaia early data release 3 (EDR3; Gaia Collaboration
et al. 2020) are cross-matched to the VPHAS+ data. Gaia EDR3 data contain
parallaxes ($\pi$), proper motions vectors in right ascension
($\mu_{\alpha}$cos$\delta$, hereafter referred to as $\mu_{\alpha}$),
declination ($\mu_{\delta}$), and the associated photometry in custom filters
measured from the first three years of Gaia observations (2014-2017) in the
2016 reference epoch. Gaia EDR3 photometry has a detection threshold of
$G\sim$ 21 mag, with the very bright ($G<$ 7 mag) and high proper-motion ($>$
0.6 arcsec yr-1) stars incomplete. The detection limit of Gaia corresponds
approximately to the 3$\sigma$ threshold of the VPHAS+ photometry in the $ri$
filters. Errors in EDR3 data release on the proper motions are believed to be
around $\pm$0.02–1.4 mas depending on source magnitude, and on parallax around
$\pm$0.02-1.3 mas. A comparison with extragalactic sources indicates an
overall negative bias of $-$17 $\mu$as (Gaia Collaboration et al., 2020) on
the parallax, which was added to our parallax measurements.

To perform cross-matching with the VPHAS+ sample, a minimum cross-match radius
of 0.1$\arcsec$ was set considering the Gaia positional accuracy
(0.059$\arcsec$), and separation between the different epochs (also
0.059$\arcsec$). Sources were initially matched using a maximum radius of
2$\arcsec$ (with both catalogues in the J2000 epoch), where the mean
separation was found to be 0.077$\arcsec$. To determine the incidence of
spurious sources, the radius was incrementally reduced by 0.05$\arcsec$, while
also counting the number of stars with magnitude differences between the Gaia
$G$ band, and the VPHAS+ $i$ magnitude (applying an offset of 0.6 mag to bring
approximately to a relation of unity) was greater than 1 mag (accounting for
variability; note that $r$-band is not used as it includes the H$\alpha$
line). A best-match radius of 0.3$\arcsec$ is chosen as 98% of matches within
1$\arcsec$ are recovered, but no sources with large magnitude offsets are
picked up. The renormalised unit weight error (RUWE) cut of RUWE$<1.4$ was
applied to remove sources with poor astrometry. As a result, 8 965 sources
with high-fidelity astrometry and photometry from Gaia EDR3 that have cross-
matches in VPHAS+ comprise our dataset. This is the dataset that will be
utilised for identifying and characterising the young stellar population
visible at optical wavelengths.

###  2.4 nIR and mIR imaging

Figure 3: ($a$): [4.5]$-$[5.8] vs. [3.6]$-$[4.5] colour-colour diagram of all
mIR sources are shown as circles. The location of sources affected by PAH
contamination and shock emission are shown following Gutermuth et al. (2009).
Stars removed following the criteria outlined in the text as contaminated
photometry are marked by filled circles. ($b$) [5.8]$-$[8.0] vs. [4.5]$-$[5.8]
colour-colour diagram of all sources are shown. Sources that are probable
external galaxies following Gutermuth et al. (2009) classification criteria
are marked by filled circles. The reddening vector for $A_{K}$=5 mag is also
shown. The mid-infrared reddening law of Indebetouw et al. (2005) was used.

The nIR sample consists of all sources with $JHK$s photometry from the UKIDSS
(UKIRT Infrared Deep Sky Survey; Lucas et al. 2008) survey. To this sample, we
discarded sources with non-stellar profiles (akin to the classification schema
in the VPHAS+ photometry), or lacking photometry in any of the three bands,
but applied no further selection criteria. This catalogue is cross-matched
with the VPHAS+–Gaia sample. Given the detection limit of UKIDSS ($\sim$19.7
mag in $J$), the optical sample misses some of the lower mass mid-M spectral
type stars (assuming average extinction, with no impact of nebulosity) that
are captured by the nIR photometry. However, nebulosity towards the core of
the Trifid nebula greatly impacts our detection of sources in the region. The
effect of this incompleteness in optical photometry compared to the nIR on our
final results is discussed in further detail in Section 4.6.

To capture the disc evolutionary stages using data from the mIR, we utilise
archival photometry from the GLIMPSE (The Galactic Legacy Infrared Midplane
Survey Extraordinaire; Churchwell et al. 2009) survey catalogue. GLIMPSE
details photometry in the [3.6], [4.5], [5.8] and [8.0] IRAC (InfraRed Array
Camera) filters along the Galactic plane, taken using the Spitzer space
telescope. Only sources with no close sources within 2$\arcsec$ and with no
photometric quality issue are considered. No photometric uncertainty limit was
applied to the data set. To this sample, we performed a cross-match search
against the MIPSGAL survey (MIPS Galactic Plane Survey; Urquhart et al. 2014)
to find any sources having photometry at 24$\mu$m, which was taken using the
MIPS (Multiband Imaging Photometer for Spitzer) instrument. Note that both the
mIR catalogues already contain cross-matches to the 2MASS $JHK$s (Cutri et
al., 2003) photometry catalogue. We cross-match the mIR catalogue with the
UKIDSS, and VPHAS+–Gaia datasets using a radius of 2$\arcsec$ (the mean
effective resolution of the mid-infrared photometry), and performed a sanity
check on the cross-matches by removing any sources which have GLIMPSE provided
$K$s band photometry deviant by more than 0.5 mag from the UKIDSS survey $K$s
photometry.

##  3 Young stellar objects

Young stellar objects (YSOs) may be discriminated from the stellar population
on the basis of their IR properties. Depending on their evolutionary stage,
YSOs are likely to be surrounded by either circumstellar discs, or envelopes
which absorb and re-radiate stellar luminosity in IR portion of the
electromagnetic spectrum. At these wavelengths, the colours of YSOs, and the
shape of their spectral energy distribution (SED) significantly differ with
respect to the vast majority of observable stars. The database of nIR and mIR
photometry spanning the wavelength range of 1.25$\mu$m–24$\mu$m described in
Section 2.4 is used to both identify YSOs, and classify their evolutionary
stage in this section.

###  3.1 Identification of young stellar objects

To identify probable YSOs, we select from the Spitzer IRAC sample only sources
with photometry in all four bands. In addition, nIR and 24$\mu$m photometry of
these sources is also used for the analysis when available. The main source of
contamination when attempting to identify YSOs against the Galactic plane is
dealing with foreground/background stars, and PolyAromatic Hydrocarbon (PAH)
and background galaxies. Most foreground/background stars do not have IR
excesses, and are centered around zero in Spitzer colour-colour planes
(Megeath et al., 2004). Removal of AGB (Asymptotic Giant Branch) stars is less
straightforward, as they do exhibit excesses similar to YSOs in some IR
colours. Similarly, background galaxies have colours falling in the region
occupied by YSOs in the mIR colour-colour plane. Active star-forming galaxies
in particular have colours that are dominated by PAH emission, and they are
bright in the [8.0] filter (due to the 7.7–8.2$\mu$m PAH feature). AGN (active
galactic nuclei), are similarly bright at 8$\mu$m and also overlap YSOs in the
colour-colour plane. In addition, the photometry also may be affected by
unresolved knots caused due to shock emission from YSOs, or from PAH emission
contaminating the apertures of sources. The latter two represent significant
sources of photometric contamination in star-forming regions. To remove extra-
galactic sources, and contaminated photometry, we adopt the criteria outlined
in Gutermuth et al. (2009), appendix A. Note that a reddening correction
corresponding to an absolute extinction, $A_{V}$ of 1.3 mag (Rho et al.,
2008), following the reddening law ($R_{V}$) of 5.5 (Cambrésy et al., 2011)
was applied. The results of this exercise are shown in Fig. 3. Overall,
$\sim$20 sources were removed following this criteria. In addition, to
identify stellar contaminants, a cutoff of [4.5]$-$[8.0]$<$0.5, and [8.0]$<5$
mag following Harvey et al. (2006) was applied. The cutoff in magnitude is
scaled to a distance of 1250 pc (see Section 4.3 for details on the chosen
distance). These sources are not removed from the sample immediately, but this
cut-off is applied when classifying YSOs.

###  3.2 Classification of young stellar objects

Figure 4: ($a$): [4.5]$-$[8.0] vs. [8.0] colour-magnitude diagram of mIR
sources are shown as circles. Sources identified as Class I and Class II YSOs
are marked as red and blue crosses respectively. ($b$) [5.8]$-$[8.0] vs.
[3.6]$-$[4.5] colour-colour diagram of all sources. The location of Class II
and Class I sources following Allen et al. (2004) and Megeath et al. (2004)
classification criteria are shown. The reddening vector for $A_{K}$=5 mag is
also shown, based on the Indebetouw et al. (2005) reddening law.

There exist multiple criteria outlined in the literature to identify and
classify YSOs using either the value of the IR SED slope ($\alpha$), or
positions of sources in colour-colour diagrams. These methods are similar in
philosophy, but differ in details. A full review and examination of these
criteria is beyond the scope of this paper, and we base the analysis in this
paper on a combination of the colour, and slope classification criteria
outlined in Megeath et al. (2004) and Allen et al. (2004), and Harvey et al.
(2006) respectively. Firstly, we classify all sources that were not identified
as extra-galactic sources, or having photometry contaminated by PAH or shock
emission using the value of $\alpha$, defined as

$\alpha=\frac{d\,log\lambda F_{\lambda}}{d\,log(\lambda)}.$ (1)

Here $\lambda$ is wavelength and $\lambda F_{\lambda}$ is flux density in
wavelength. To determine $\alpha$, we fit the fluxes of all sources longward
of 2$\mu$m. For all sources, we require at least photometry in all IRAC bands.
The photometry was dereddened assuming the mean extinction towards the region.
The slope was calculated using a weighted least squares fit to the available
photometry. Here, the weight applied was the inverse of the squared
photometric errors. The fit error from a standard weighted least squares fit
is also calculated.

Following the classification scheme of Harvey et al. (2006), we classify
sources as YSOs based on their $\alpha$ value. The YSO classes are;

  1. Class I: $\alpha>$ 0.3

  2. Flat: $0.3\geq\alpha\geq-$0.3

  3. Class II: $-0.3>\alpha\geq-$1.6

Note that we do not attempt to classify Class III sources. Class III sources
(YSOs with no or transition discs) cannot be differentiated from the bulk of
the stellar population as they lack IR excesses. While in the literature Class
III sources have been classified as non-contaminants with $\alpha<-$1.6, we
refrain from doing so here as the Trifid nebula lies in a direct sight line
towards the Galactic plane, and even in mIR is heavily contaminated by
background/foreground field stars (Feigelson et al., 2013). Therefore,
classifying Class III sources based on $\alpha$ would lead to spurious
detections.

Spitzer images of the sources were visually inspected to detect any that may
be affected by extended emission, nearby contaminating bright sources (recall
that we discard previously any sources having counterparts within 2$\arcsec$,
the mean full-width half maximum of the photometry), or are elongated. Three
sources that are extended, or have a nearby bright source within 2$\arcsec$5
are removed. Note that the bulk of our sample has $\alpha\sim-2.5$ to $-3$,
agreeing with the expectation for main-sequence stars.

Table 1: Infrared photometry of young stellar objects in the Trifid Nebula.

GLIMPSE ID | $J$ | $H$ |  $K$s | [3.6] | [4.5] | [5.8] | [8.0] | $\alpha$ | Class
---|---|---|---|---|---|---|---|---|---
| (mag) | (mag) | (mag) | (mag) | (mag) | (mag) | (mag) |  |
G006.9115-00.2629 |  |  |  | 12.07$\pm$0.06 | 10.30$\pm$0.04 | 9.26$\pm$0.04 | 9.18$\pm$0.09 | 1.02$\pm$0.3 | I
G006.9258-00.2629 | 13.27$\pm$0.02 | 12.33$\pm$0.02 | 11.73$\pm$0.033 | 10.88$\pm$0.05 | 10.55$\pm$0.04 | 10.17$\pm$0.06 | 9.08$\pm$0.06 |  $-$0.9$\pm$0.2 | II
G006.9292-00.2460 | 16.30$\pm$0.10 | 14.08$\pm$0.05 | 12.13$\pm$0.02 | 9.53$\pm$0.03 | 8.46$\pm$0.03 | 7.60$\pm$0.03 | 6.76$\pm$0.02 | 0.2$\pm$0.1 | Flat
G006.9451-00.3318 |  | 14.18$\pm$0.07 | 12.45$\pm$0.05 | 10.60$\pm$0.17 | 10.00$\pm$0.13 | 9.32$\pm$0.04 | 8.96$\pm$0.03 |  $-$1.49$\pm$0.3 | II
G006.9637-00.3104 | 11.15$\pm$0.02 | 10.75$\pm$0.02 | 10.42$\pm$0.02 | 9.7$\pm$0.08 | 9.23$\pm$0.04 | 8.26$\pm$0.082 | 6.97$\pm$0.18 | 0.27$\pm$0.4 | Flat

Note. — Five selected rows are shown to display the form and contents of the
table. The full table is only available electronically, along with the source
coordinates. The 24$\mu$m magnitude for G007.0317-00.2845 is 1.7$\pm$0.03,
G006.9221-00.2513 is 2.29$\pm$0.02, G007.0097-00.2542 is 2.87$\pm$0.18,
G006.9292-00.2460 is 4.14$\pm$0.03, G007.0574-00.2701 is 4.79$\pm$0.12 mag
respectively.

After the SED slope analysis, an additional colour selection is applied to cut
any remaining contaminants (e.g., AGB stars, reddened main-sequence stars).
The criteria outlined in Megeath et al. (2004) and Allen et al. (2004) is
utilised, and the results are shown in Fig. 4. Here, the approximate colours
of Class II, Class I, a combination of Class I/II sources, and reddened Class
II YSOs given by Megeath et al. (2004) and Allen et al. (2004) are shown. The
location of sources classified based on $\alpha$ are also shown. The loci of
main sequence sources without excess in the IR are centred around zero as
expected, with a spread due to reddening differences and CO absorbing giants.
The positions of Class II sources are demarcated by the bounding box, which is
reproduced by models of discs with varying accretion rates, and inclinations
around young and low-mass stars (D’Alessio et al., 2001). The empirical
boundary between Class II and Class I sources is also given, and their colours
are of protoestellar objects with infalling envelopes. The location of Class
II sources affected by [8.0] excess emission is marked in the region of Class
I/Class II sources, as is the location of Class II sources affected by extreme
reddening. While in general, the location of these sources agrees with the
expectation based on their colour criteria, we note that 17 sources fall
blueward of the region demarcated as either Class II or reddened Class II
sources. Sources falling in this location are spatially preferentially located
towards the edges of the nebula, directly in the sight of the Galactic plane,
and have colours [4.5]$-$[8.0]$<$1 (Harvey et al., 2006), similar to AGB
stars. These are likely reddened AGB stars in the Galactic plane. To remove
them, these 17 sources falling outside of the expected colour boundaries for
YSOs from Megeath et al. (2004) and Allen et al. (2004) are removed.

###  3.3 Final sample of YSOs

Our final sample consists of 46 Class II YSOs, 9 Flat, and 6 Class I YSOs. The
complete sample along with archival photometry is given in Table 1. A source
of contamination in our final sample is background/foreground AGB stars, as
sources resembling colours of galaxies, or contaminated aperture photometry
were removed, and the images of the final classified sources were visually
inspected. The contamination rate in our sample may be judged by adopting the
locations of AGB stars in colour-magnitude diagrams as found by Harvey et al.
(2006). In that paper, it is suggested that some AGB stars which may have
similar SED slopes to YSOs have colours [4.5]$-$[8.0]$<1$. These AGB stars
represent the largest source of contamination in our sample, given the sight
line towards the Galaxy. 12 YSOs have colours bluer that AGB stars (but redder
than [4.5]$-$[8.0]$>$0.5), and these represent $\sim$20% of our sample. This
represents an upper limit on the estimate of the contamination in the sample,
given that otherwise the sources meet the SED and colour classification
criteria.

The completeness of the sample is difficult to judge straightforwardly, as
there are various factors affecting it. In particular, the saturation of
photometry towards the centre of the Trifid, and the large-scale differences
in the nebulosity. These cannot be straightforwardly estimated even using a
control field given the location near the Galactic bulge. A similar conclusion
was arrived at by the study of Feigelson et al. (2013), who identified YSOs in
the Trifid nebula based on a combination of modelling IR SEDs, and X-ray
emission.

As an alternative, we estimate the recovery rate of our sample against the one
selected by Feigelson et al. (2013). They identified 41 Stage II/III, and 22
Stage I YSOs having Spitzer photometry in all IRAC bands within the spatial
boundary used in this work by modelling the available photometry against the
YSO models of Robitaille et al. (2006). From their sample, we find 28 Stage
II/III, and 11 Stage I YSOs have cross matches in our final YSO sample. In
their remaining Stage II/III sample, we find that two have mIR colours
resembling PAH contaminated sources, while two are likely stellar contaminants
([4.5]-[8.0]$<$0.5) according to our analysis, and three had $\alpha$ values
of Class II YSOs, but were removed from the final sample based on their
positions in the colour-colour plane resembling contaminants (see Section
3.2). The remaining six had $\alpha$ values between $-1.7$ and $-2.1$ (two of
which fell in the contaminant region in the colour-colour plane). This
distribution is expected given that the Stage and SED slope classifications
are not analogous, and depend on various factors (Robitaille et al., 2006). In
addition, the Stage classification assumes that sources with $\alpha>-2$ are
Stage II, and $\alpha<-2$ are Stage III. Out of the 11 Stage I sample from
Feigelson et al. (2013) not having counterparts in this sample, 1 star was not
present in our source catalogue, while 7 were discarded as PAH or shock
emission contaminants. The three remaining sources have colours resembling
reddened AGB stars and were removed from our final sample. Therefore, our
sample compares well to the YSO catalogue of Feigelson et al. (2013) once the
constraints chosen in this work are applied.

Given the contamination rate, and absence of a $\chi^{2}$ fitting value (to
adjudicate the quality of the classification) for the sources in Feigelson et
al. (2013), archival infrared Class III sources from that catalogue are not
included here. Similarly, X-ray bright Class III sources from that survey were
not included as the X-ray detection threshold chosen is considerably faint,
and a significant fraction ($\sim$35%) of Class III X-ray sources do not have
counterparts at other wavelengths, suggesting they could be faint background
sources. In addition, a distance value of 2700 pc was (see Section 4.3) used
to calculate membership probabilities and absorption, possibly leading to the
inclusion of more distant X-ray sources.

##  4 Pre-main sequence stars

As YSOs evolve, their surrounding material collapses into circumstellar discs
to conserve angular momentum. Accretion of mass from the disc towards the
central stellar source via the stellar magnetosphere is necessary to gain
enough material to initiate fusion. This accretion process results in a
combination of unique signatures including excess ultraviolet continuum, and
line emission, notably in H$\alpha$ (Gullbring et al., 1998). It is through
this excess H$\alpha$ emission that accreting young PMS stars can be readily
identified.

Figure 5: ($r-i$) vs. ($r-$H$\alpha$) colour-colour diagram of all sources
considered in the work, shown as dots. A reddening vector of 1 mag is shown in
the top left. The dashed line represents the main-sequence colour locus
calculated from the Pickles (1998) synthetic spectra, reddened by $A_{V}=1.3$
mag. Stars selected as accretors based on the EWHα criterion are marked by
blue crosses. Figure 6: ($r-i$)0 vs. EWHα diagram of all sources. Stars
selected as accretors are shown as blue crosses, along with their associated
errors. The selection criteria is shown as a dashed line, and the spectral
type at different values of ($r-i$)0 are marked.

###  4.1 Identification of accreting pre-main sequence stars based on
photometric EW${}_{\textrm{H}\alpha}$

It has been established that the ($r-i$) vs. ($r-$H$\alpha$) diagram (Fig. 5)
can be used to obtain estimates of the H$\alpha$ equivalent width (EWHα) of
stars with precision photometry (De Marchi et al., 2010), given the extinction
and spectral type are relatively well-approximated (to within 0.2 mag). We
follow the method detailed in Kalari (2019) for VPHAS+ photometry to estimate
the H$\alpha$ equivalent width of all VPHAS+ sources, adopting the average
absolute visual extinction value ($A_{V}$) of 1.3 mag from Rho et al. (2008).

From the values of photometrically estimated EWHα, we select accreting PMS
candidate stars based on their ($r-i$) colour (as a spectral type proxy). Note
that ($r-i$) colours, and $r$-band magnitudes are corrected for H$\alpha$
emission, as the H$\alpha$ line falls within the $r$-band filter. Finally,
images of accretors were also visually inspected using unsharp masking to
remove sources affected by nebulosity. The reader is referred to Kalari et al.
(2015) for details on these procedures, and also the H$\alpha$ excess
candidate selection criteria. To select PMS candidates, we use the spectral
type-EWHα criteria for Classical T-Tauri stars (CTTS) of Barrado y Navascués &
Martín (2003), adding to the upper limit of selected candidates based on the
average error values of EWHα (CTTS are low-mass accreting PMS stars, showing
H$\alpha$ in emission). Stars earlier than K5 spectral type having EWHα values
less than $-12$Å (where negative values denote in emission), in the K5-M2
spectral range with EWHα values less than $-$15Å , and stars later than M2
with H$\alpha$ in emission less than $-$25Å are selected as candidate PMS
stars. Results of the selection procedure can be seen in Fig.6. The selection
criteria account for propagated errors from the photometry, and uncertainty in
spectral type. They are also large enough to weed out main-sequence
contaminants with solely chromospheric H$\alpha$ emission. 68 PMS candidates
were selected based on this method. The resulting EWHα values for the selected
candidates, along with their photometry and astrometry is given in Table 2.

###  4.2 Kinematic properties of pre-main sequence candidates

Gaia proper motions provide a valuable membership indicator for nearby star-
forming regions. Given that the Trifid nebula is located in our line of sight
towards the plane, it is likely that non-members may be erroneously identified
as PMS members. Using proper motions, these can be excluded from further
study.

M 20 is identified in the $\mu_{\alpha}$–$\mu_{\delta}$ plane based on the
mean motion of PMS candidates, having errors in any proper motion vector less
than 0.5 mas yr-1. The cluster itself is separated in proper motion space with
respect to the surrounding field population, with an increase in density
visible around the cluster centre identified based on the motions of the PMS
candidates (see Fig. 7), and verified through the density map. Notice the
large scatter of sources around the mean proper motion. This is expected given
the location of the region in the Milky Way, as the photometry contains
numerous background/foreground stars. The exact centre of the region can be
best distinguished using the identified PMS candidates.

To calculate the centroid in the proper motion plane, we fit double Gaussians
to the proper motion vector histograms, which represent a conflation of the
field population, having a broad distribution, and the cluster population,
which has a well-defined peak. The proper motion centres in both directions
was estimated based on the narrow Gaussian parameters. For $\mu_{\alpha}$, the
centre is 0.22$\pm$0.8 mas yr-1, and in $\mu_{\delta}$ it is at
$-$1.71$\pm$1.6 mas yr-1. The error bar denotes the width of the Gaussian
peak.

Proper motion candidates are selected based on the estimated proper motion
centre of the cluster. We assume a median width of 2 mas yr-1 from the cluster
centre. PMS candidates meeting the proper motion criteria are selected as
those lying between $-1.8<\mu_{\alpha}<2.2$ and $-3.7<\mu_{\delta}<0.3$. Stars
lying outside this bounding box are rejected as members, if their proper
motions values along with the errors lie outside this boundary. 27 candidates
were removed based on this criterion. Since the proper motion centre of the
PMS members lies within the field proper motion distribution, it is possible
that some non-members may have been included in our final analysis.

Finally, non-accreting members that can potentially be identified using Gaia
proper motions are not included as the proper motion distribution of the
cluster overlaps significantly with the field. A selection criterion to
isolate such stars using only proper motions would include a non-negligible
fraction of non-members, with no further defining physical characteristic to
identify members.

Figure 7: Gaia EDR3 proper motions of PMS candidates in the Trifid Nebula
(green circles). The underlying inverted grayscale density map reflects the
distribution of all Gaia sources in the field of view. Candidates selected as
members on the basis of their proper motions are shown by blue crosses, while
those rejected as members based on the proper motion criteria are shown as red
circles. The proper motion selection criteria is represented by the solid
bounding box. Dots represent proper motions of all stars in the area of study.

###  4.3 Distance to the Trifid nebula

A large discrepancy exists in the literature on the distance towards M 20.
Previous distance measurements have been made using spectroscopic parallax
(Rho et al., 2008; Tapia et al., 2018), and extinction measurements (Cambrésy
et al., 2011). Spectroscopic parallax measurements of HD 164492A place it at a
distance between $\sim$1.6–2.0 kpc (Rho et al., 2008; Tapia et al., 2018),
depending on the reddening law chosen (see Tapia et al. 2018). This places it
behind the Sagittarius arm, and in front of the Scutum arm of the Galaxy.
Whereas distance measurements using extinction mapping by Cambrésy et al.
(2011) suggest a much farther distance of around 2.7 kpc, placing it firmly in
the Scutum arm of the Galaxy. In this paper, we eschew from depending on a
priori knowledge of the reddening towards the region to measure the distances,
by using high fidelity parallax values of identified PMS candidates.

To estimate the distance towards the region, we consider only PMS candidates
with Gaia EDR3 $\pi$ errors ($\sigma_{\pi}$) less than 0.1 mas. As this
criterion leaves us with only three stars, we add to the sample five
previously identified members from the YSO sample, which all have
$\sigma_{\pi}<$ 0.1 mas. This increase from the IR identified YSO sample is
expected as that sample is concentrated towards brighter magnitudes, where the
$\sigma_{\pi}$ are lower. In total, we have 8 stars that are candidate members
of the region, with high-fidelity $\pi$. Our chosen error criteria is to
remove faint stars with non-negligible errors, that may degrade the accuracy
and precision of the final distance estimate.

To estimate the individual distances of the chosen sample, we used the method
of Bailer-Jones et al. (2018). An inference approach is adopted, which uses a
prior with a length scale parameter (chosen here to be 1750 pc, based on the
Galactic model given in Bailer-Jones et al. 2018). The resulting confidence
intervals of the individual distances are asymmetric. Based on this method the
resulting distance histogram of all candidates is shown in Fig. 8. An inverse
weight of the mean error was applied to the histogram. A Gaussian is fit to
the histogram, accounting for the well defined cluster distance. From the peak
of the Gaussian fit, we estimate a distance of 1257${}^{+190}_{-98}$ pc to M
20 (the median distance of the sample is 1255 pc). The distance towards the
cluster is taken to be 1250 pc for our analysis. The asymmetric error
represents the median of the 5th and 95th percentile confidence intervals of
the resultant distance distribution.

The resulting distance from Gaia EDR3 parallaxes is much closer than previous
measurements using spectroscopic parallaxes or extinction mapping. It agrees
with the estimate of $\sim$1.22 kpc made using Gaia DR2 parallaxes by Wade et
al. (2017). As a comparison, also shown are the histogram of distances of Gaia
EDR3 sources meeting the proper motion criterion, and having RUWE and
$\sigma_{\pi}$ within the limit for our sample. This sample has a median
distance of 1300 pc, but may contain few contaminant background stars.
Finally, the new closer distance places the cluster firmly in the Sagittarius
arm, alongside the star forming regions in the near vicinity such as the
Lagoon Nebula, Simeis 188 (both lying below the Galactic plane), and the
optically invisible H II region W28 A2, and its supernova remnant W28.

Figure 8: Histogram of distances of candidate members having
$\sigma_{\pi}<$0.1 mas (hatched bars). The resulting Gaussian is shown as a
dashed line, and has a peak of 1250 pc. Also shown only for comparison are the
distances of all Gaia EDR3 sources meeting the proper motion membership and
parallax error criteria.

###  4.4 Stellar properties of pre-main sequence candidates

We estimate the stellar properties (mass, $M_{\ast}$ and age, $t_{\ast}$) of
the 41 PMS candidates that meet the proper motion criteria, by interpolating
their positions in the observed ($r-i$) vs. $r$ colour-magnitude diagram
compared to model isochrones and mass tracks (Fig. 9). The stellar isochrones
and mass tracks of Bressan et al. (2012) in the VPHAS+ filter set (AB
magnitude system) are utilised for this purpose. The AB to Vega offset values
for the $r$ and $i$-band filters were 0.174 and 0.378 respectively. To
understand the differences between the choice of stellar models, the
$M_{\ast}$ and $t_{\ast}$ are compared to those estimated using the Siess et
al. (2000) and Dotter et al. (2008) models in the VPHAS+ filter set. The
metallicity in the models is set to solar, and the models only account for
single stars. The $A_{V}$ is assumed to be uniform, and is set at 1.3 mag,
with an $R_{V}$ of 5.5 (Cambrésy et al., 2011) used to deredden the observed
photometry. A distance of 1250 pc was used to transform the models to the
observed distance. The results are reported in Table 2. Errors on the
interpolated stellar properties include propagated photometric uncertainties,
and an assumed distance and extinction uncertainty of 200 pc, and 0.2 mag
respectively. The effect of assuming a uniform extinction value on our results
is also discussed.

Table 2: Optical photometry, Gaia astrometry, and estimated stellar properties
of pre-main sequence stars in the Trifid Nebula.

VPHAS+ ID | $r$ | $i$ | H$\alpha$ | EWHα | $\mu_{\alpha}$ | $\mu_{\delta}$ | Mass | Age
---|---|---|---|---|---|---|---|---
| (mag) | (mag) | (mag) | (Å) | (mas yr-1) | (mas yr-1) | ($M_{\odot}$) | (Myr)
816177248012 | 20.07$\pm$0.07 | 18.53$\pm$0.03 | 18.53$\pm$0.06 |  $-$149.42$\pm$11.1 |  $-$1.173$\pm$0.678 |  $-$1.508$\pm$0.412 | 0.19${}^{+0.01}_{-0.02}$ | 2.24${}^{+0.39}_{-0.44}$
816177242912 | 13.41$\pm$0.01 | 12.99$\pm$0.01 | 13.11$\pm$0.01 |  $-$14.56$\pm$1.5 |  $-$0.011$\pm$0.596 |  $-$3.23$\pm$0.433 | 1.88${}^{+0.36}_{-0.23}$ | 6.08${}^{+2.73}_{-3.00}$
816177256886 | 20.19$\pm$0.07 | 18.47$\pm$0.03 | 18.96$\pm$0.07 |  $-$70.25$\pm$4.1 | 0.522$\pm$0.544 |  $-$0.522$\pm$0.342 | 0.18${}^{+0.01}_{-0.02}$ | 1.81${}^{+0.28}_{-0.35}$
816856488067 | 19.55$\pm$0.044 | 18.25$\pm$0.025 | 18.23$\pm$0.035 |  $-$113.88$\pm$6.5 |  $-$0.667$\pm$0.538 |  $-$2.75$\pm$0.319 | 0.25${}^{+0.01}_{-0.02}$ | 3.01${}^{+0.68}_{-0.8}$
816856492037 | 19.87$\pm$0.05 | 18.18$\pm$0.02 | 18.79$\pm$0.05 |  $-$46.22$\pm$2.6 | 0.123$\pm$0.529 |  $-$0.609$\pm$0.313 | 0.21${}^{+0.01}_{-0.01}$ | 1.58${}^{+0.26}_{-0.37}$

Note. — Five selected rows are shown to display the form and contents of the
table. The full table is only available electronically, along with the source
coordinates. Mass and age errors denote the upper and lower deviations based
on the propagated photometric uncertainties, and assuming a distance and
reddening uncertainty of 200 pc, and 0.2 mag.

The resulting colour-magnitude diagram of the PMS stars, with the Bressan et
al. (2012) stellar tracks and isochrones overlaid is shown in Fig. 9. PMS
candidates from our sample fall approximately between the 0.2–2.0 $M_{\odot}$
tracks, and are clustered between the 1 and 2 Myr isochrones. The results
estimated from the three different stellar models are given in Fig. 10. Here,
the median mass of our candidates peaks around 0.3 $M_{\odot}$, with the most
massive stars $\sim$ 2 $M_{\odot}$. This mass distribution is expected given
the absolute $r$-band luminosity function of our candidates (Fig. 10c), where
our candidates fall between 2–9 mag (after correcting for extinction and
distance). This approximates to early F– early M candidates, which corresponds
to the derived mass range. Also marked for reference in Fig. 10c are the
expected absolute $r$-band magnitudes for main-sequence stars. There are no
significant differences between stellar masses due to the assumption of a
particular set of models accounting for standard errors. Note that the mass of
the most massive stars according to the Siess et al. (2000) models is slightly
higher when compared to the other two model results.

The median age of the PMS candidates interpolated using the Bressan et al.
(2012) isochrones is 1.5${}^{+0.3}_{-0.5}$ Myr. Error bars correspond to the
median values for the upper and lower age errors, which are estimated from
interpolated ages including the photometric, reddening, and distance
uncertainties. The age distribution of the PMS candidates is shown in Fig.
10a. Note that the reddening vector exhibits a similar slope to the isochrones
in the ($r-i$) vs. $r$ plane (Fig. 9). Minor increments in extinction will not
change significantly the estimated age of the region, as evidenced by the
resulting age errors. The median ages from the Bressan et al. (2012), Siess et
al. (2000), and Dotter et al. (2008) models are 1.5, 1.4, and 1.9 Myr
respectively, and that the estimated age differences between models are lower
than the standard errors. For the remainder of this paper, we utilise the
stellar properties derived using the Bressan et al. (2012) models.

The median age of the PMS stars estimated ($\sim$1.5 Myr) here is older than
the literature age of the region $\sim$0.3 Myr. This age estimate was made
based on the expected wind-bubble shell size by Cernicharo et al. (1998).
Although this value does not take into account earlier stages, it represents a
good lower limit to the cluster age. The photo-ionisation time of cometary
globules in the region (Lefloch et al., 2002), and the presence of dense star-
forming cores (Tapia et al., 2018) attest to a young age ($\lesssim$1 Myr), as
being the most recent burst of star formation in the region (Rho et al.,
2008). The age estimated in this work represents well only the current
generation of PMS stars in the region. This generation of stars must have
formed alongside or slightly before HD 164492A, as its spectroscopic age is
around 0.6 Myr (Petit et al., 2019).

The impact of our assumption of a uniform extinction value $A_{V}$=1.3 mag on
our resultant ages is tested by adopting higher uniform values of extinction.
Assuming $A_{V}$=2 mag, increases the median age of the sample slightly to
$\sim$1.8 Myr. This agrees with the expectation of the reddening vector
discussed earlier. Assuming a higher extinction of $A_{V}$=5 mag, around 10
stars fall beyond the 100 Myr isochrone, and their ages cannot be estimated.
These stars would be considerably blue, and lie beyond the main-sequence at
the distance of 1250 pc. The median age of the remaining sample is 2.8 Myr.
Thus, it is likely that the PMS stars have a lower mean extinction than
assumed here, than higher. Finally, we found that a simulation of the random
variation of individual stellar reddening for PMS stars leads to an age
estimate of 1.7 Myr. To arrive at this value, we generated random values for
the assumed extinction correction (between $A_{V}$=0.5 and 5 mag) for each
star. Each star was also cloned 10000 times, and the median age of these
clones with random extinction correction was the result. These experiments
serve to demonstrate that the assumed mean extinction is not a significant
under-estimation for the majority of our sample; and that while individual
stellar extinction values are essential for a precise age estimate, the final
median age of our sample should not vary significantly. From this analysis, we
also suggest that future spectroscopic surveys of the young stellar population
that can estimate precise individual extinctions and spectroscopic properties
can pin down the precise age of individual stars in the Trifid nebula.

Figure 9: ($r-i$) vs. $r$ colour-magnitude diagram of PMS candidates.
Overplotted and labelled are the isochrones, and mass tracks from (Bressan et
al., 2012). The extinction vector for $A_{V}=1$ mag is shown in the top right.
Figure 10: Age, Mass, and $M_{r}$-band histograms of PMS candidates in the
Trifid Nebula in panels (a), (b), and (c) respectively. For the age, and mass
histograms, the solid (blue), dotted (green), and dashed (red) lines represent
values estimated using the Bressan et al. (2012), Dotter et al. (2008), and
Siess et al. (2000) models respectively. In the $M_{r}$ histogram the
approximate absolute $r$-band magnitude of select spectral types from Pecaut &
Mamajek (2013) are labelled.  Figure 11: (a): ($H-K$s) vs. ($J-H$) colour-
colour diagram of the PMS sample. The reddening vector for $A_{V}$ =1 mag is
also shown in the top left. The CTTS locus from Meyer et al. (1997) is shown
as a dashed line. The main-sequence colours from Pecaut & Mamajek (2013) are
shown by a dotted line. (b): The mIR colours of the sample are shown. The
boundary of colours for Class II sources with a median accretion of 10-8
$M_{\odot}\,yr^{-1}$ from the models of D’Alessio et al. (2001) are shown,
along with the positions of Class I and Class III sources. In both panels, the
the size of the symbol is scaled with the estimated stellar age, given in the
index in the lower right of panel (a).

###  4.5 Infrared properties of pre-main sequence candidates

As described in Section 3, circumstellar discs surrounding PMS stars are
heated by the central star. This leads to an inner disc wall which re-radiates
and absorbs energy primarily at nIR wavelengths leading to prominent $K$-band
excesses compared to main-sequence stars. Evidence for nIR–mIR excesses in PMS
stars are thus excellent independent indicators of circumstellar discs, and
accretion. The exact value of the IR excesses depends on many factors (Meyer
et al., 1997). Note that an absence of excess does not conclusively rule out
the presence of a cirumstellar disc, but could be due to an inner disc hole,
or edge-on disc inclination angle.

To identify if our candidate PMS stars have near or mIR excesses resembling
circumstellar discs, we use the cross-matched nIR and mIR data described in
Section 2.4. 39 PMS candidates have nIR counterparts in $JHK$s bands from the
UKIDSS survey meeting the photometric criteria. The ($H-K$s) vs. ($J-H$)
colour-colour diagram is shown in Fig 11a. The main-sequence locus from Pecaut
& Mamajek (2013) is shown along with the extinction vector for 1 mag. Also
marked is the CTTS locus from Meyer et al. (1997), showing the expected nIR
excesses for CTTS stars having accretion rates between 106–108
$M_{\odot}\,yr^{-1}$. 90 % of our sample have colours within errors falling on
the CTTS locus, suggestive that the majority of our sample have circumstellar
discs. The lack of nIR excesses for a few stars in our sample does not suggest
the absence of circumstellar discs, and could be due to inclination angle or
inner disc hole.

While 32 PMS candidates have counterparts in one or more mIR Spitzer band,
only eight have photometry in all bands and are shown in Fig. 11b. Seven are
marked as Class II YSOs from the $\alpha$ analysis (with the remaining star
having $\alpha=-1.9$). For the GLIMPSE survey, the limiting magnitude in the
Spitzer [3.6] band is around 14 mag, which is approximately a mid-K spectral
type. Therefore, it is not surprising that while nearly all our candidates
have nIR counterparts with high quality photometry from UKIDSS, only a
brighter sub-sample have mIR photometry. The expected colours for stars having
mean accretion rates of 10-8 $M_{\odot}\,yr^{-1}$ from D’Alessio et al. (2001)
is shown, along with the expected positions of Class I and Class III sources.
All of the sample falls in the expected positions for Class II sources with
the mean accretion rate around 10-8 $M_{\odot}\,yr^{-1}$. The nIR and mIR
colours of our sample clearly suggest that they are PMS stars. The IR colour-
colour diagrams shown in Fig. 11 also relates the position in colour space to
the estimated age. No clearly demarcated trend is seen from the observed
colours. Stars that are detected at mIR wavelengths are younger than the
general sample. The photometry of the PMS sample in nIR and mIR is given in
Table 2.

###  4.6 Sample completeness

Sample completeness is a function of photometry meeting our quality criterion
for all sources in the Trifid nebula in magnitude and spatial location.
Completeness in magnitude space estimated from the luminosity function in
$ri$H$\alpha$ are 19.8, 18.5, and 19.4 mag respectively (see Kalari et al.
2015 for details). Based on the photometric magnitudes, we are complete until
a stellar mass around 0.45 $M_{\odot}$. Spatial variation of completeness is
considerable in our sample due to nebulosity varying on small spatial scales.
In this scenario, high-fidelity photometry of sources cannot be determined
using currently available sky subtraction techniques (see Section 2.2). Note
that while these sources are identified, their photometry particularly in
H$\alpha$ is sufficiently degraded to be unsuitable for further analysis. To
gauge this incompleteness, we estimate the fraction of the total number of
sources detected in the $ri$H$\alpha$ magnitude across the spatial extent of
the study, to the number of sources meeting the photometric quality criteria
in equally spaced grids. This fraction is represented in Fig. 12. This is
created by assuming that for an increasing aperture radius, if the increase in
observed magnitude is larger than the error on the magnitudes, the photometry
is affected by the increasing nebulosity (i.e., being classed as an extended
source in Drew et al. 2014). This is given by the equation

$|r_{\textrm{4}}-r_{\textrm{3}}|>3\times\sqrt{{r_{\textrm{4}}}^{2}+{r_{\textrm{3}}}^{2}}+0.05.$
(2)

Here $r_{\textrm{4}}$ and $r_{\textrm{3}}$ represent the magnitudes in
apertures of 2.8$\arcsec$ and 2$\arcsec$ respectively. The smoothed density
map of sources affected by nebular subtraction in Fig. 12 suggests that the
sample is more incomplete towards the central dust lane, and HD 164492A.
Accounting for the missing fraction of PMS stars due to this spatial
incompleteness is difficult, as it depends on the inherent spatial
distribution, and their fraction compared to the total number of sources in a
given region. This caveat is considered when analysing the spatial
distribution in the Trifid in the region around HD 164492A. Also over plotted
in Fig. 12 is the spatial distribution of Gaia EDR3 point sources, as a
comparative reference for high-quality optical imaging in the region.

Overall, we can conclude that our sample is complete down to masses of 0.45
$M_{\odot}$, and spatially our sample is significantly affected by the
nebulosity. The effect of this incompleteness suggests sources missing around
the central dust lane, and near HD 164492A. Based on this, we can assume that
our sample represents only a lower limit of PMS stars around HD 164492A. Given
these circumstances, we define a radius around the central O star where the
sample is incomplete, and consider this caveat when discussing the spatial
properties of the PMS sample.

Figure 12: Colour map in right ascension and declination depicting ratio of
sources detected to those whose photometry is affected by spatially varying
nebulosity on small scales. Gaia EDR3 sources are overlaid as grey circles.

##  5 Discussion

###  5.1 Star formation triggered by CCC

As discussed in Section 1, the Trifid Nebula is a proposed site of CCC
triggering star formation. From a sub-mm study of CO lines Torii et al. (2011,
2017) identified three distinct velocity components at 1.4, 9 (here after 2 km
s-1 cloud, and cloud C respectively following their nomenclature), and 18.2 km
s-1, representing three different molecular clouds. They claimed that the 2 km
s-1 cloud, and cloud C are observed to have higher temperatures, and are
likely heated by the central ionising source (HD 164492A) based on positional
and luminosity arguments– suggesting that they are the parent clouds of the
central cluster (while the 18 km s-1 cloud is unlikely to be related to
current star formation in the region). Torii et al. (2011) indicate that the
observed velocity pattern is systemic. The 2 km s-1 cloud is presently moving
towards us, with cloud C moving away from us with an observed velocity
separation of 7.5 km s-1. Here, the velocity separation refers to the
differences in the peak velocity of the two clouds CO spectra, and is taken
from Torii et al. (2011). The 2 km s-1 cloud is associated with the dark lanes
in the Trifid nebula, but cloud C does not correspond to the optically
observed features strengthening the argument that cloud C lies behind the 2km
s-1 cloud in our line of sight. The total estimated molecular and stellar mass
is insufficient to gravitationally bind the system, leading the authors to
propose a different mode of star formation. They suggest that a consistent
scenario with the observed velocity and spatial properties of the molecular
clouds is that these two clouds collided $\sim$1 Myr ago, with the observed
relative velocity difference arising from the collision (see also Fig.13). The
collision between these two clouds triggered rapid star formation at the
location of the collision, including the formation of HD 164492A.

Important differences between stars forming in a single molecular cloud versus
a pair of colliding clouds are the spatial positions and age spreads of the
newly formed stars with respect to the position and velocities of their natal
molecular clouds (Loren, 1976). When gravitational collapse of molecular
clouds is induced by turbulence, it is thought that heavier filaments fall
towards the centre, and the majority of stars in the central region form over
a prolonged period of time from such dense material giving rise to a
concentration of young stars in cluster centres, with the more massive stars
forming there preferentially. Older stars are expected to be detected more
towards the outskirts of regions either due to lower molecular densities in
the outer areas, resulting in the earliest period of star formation exhausting
material to form stars, or due to dynamical interactions over time. In the
case of head-on CCC as considered in M 20, the new generation of stars will
form where the densities and compression are the highest. This is expected to
be towards the edges of the collided clouds, as a cavity is created towards
the centre of the collision, along with a detected bridge feature at
intermediate velocities (Torii et al., 2017). The stars which formed rapidly
after the collision will now be found in a region where the two separate
velocity components are coalescing into one intermediate velocity. The younger
generation of star formation currently occurring will be concentrated towards
the edges of the collided clouds, and near the edges of the cavity where dense
molecular material is not yet exhausted. Stars formed as a result of the
collision will be found exclusively around the collapsing region, but the
central overlapping region where the clouds collided might be devoid of stars.

Therefore, the spatial and age spreads of young stars can be correlated with
the velocity and spatial properties of molecular clouds and filaments to
identify the sites of CCC induced star formation. Based on the demographics of
the young stellar population identified in Sections 3 and 4, we explore here
whether they provide evidence for or against the CCC hypothesis.

Figure 13: Herschel 350$\mu$m/8$\mu$m Spitzer/VPHAS+ H$\alpha$ $rgb$ image of
the Trifid Nebula. North is up and east is left. The scalebar corresponding to
1 pc at the distance of the Trifid nebula is also given. Circles mark the
location of PMS stars detected in our sample. Their size is related to their
age, and scaled according to the legend in the top left corner. Blue and cyan
crosses mark the position of Class II and Class I sources respectively, while
the plus sign indicates the location of star-forming cores from Tapia et al.
(2018). The blue and red contours indicate the blue-shifted 2 km s-1 cloud,
and red-shifted cloud C which have collided following Torii et al. (2011). The
contour levels are drawn for clarity at 18 and 36 K km s-1 for both clouds.
Black contours represent the filaments identified from Herschel imaging.

####  5.1.1 State of colliding clouds in Trifid nebula

The mass of the two clouds is calculated based on the CO luminosities and line
widths reported in Torii et al. (2011), where the authors adopted the $X_{\rm
CO}$ factor of 2.0 $\times$ 1020 cm-2(K km s-1)-1 to convert the CO luminosity
to mass. The difference between the calculations reported here and those in
Torii et al. (2011) is in the distance adopted. For the calculations reported
in this paper the distance adopted is 1250 pc from Gaia EDR3 parallaxes,
wheras Torii et al. (2011) adopted much larger distance estimates of 1700 and
2700 pc. Based on the new distance, the total mass and radius of cloud C is
estimated as 0.5$\times 10^{3}$ $M_{\odot}$ and 0.9 pc respectively. The 2 km
s-1 cloud has a mass of 0.35$\times 10^{3}$ $M_{\odot}$ and a radius of 0.7
pc. Considering the gravitational situation of the cloud, we can calculate
whether it is bound or unbound based on the $\alpha_{\rm G}$ parameter. If
$\alpha_{\rm G}>1$ the cloud is considered unbound, and if $\alpha_{\rm G}<1$
the cloud is bound. Here,

$\alpha_{\rm G}=\frac{5\sigma_{v}\,R}{GM}.$ (3)

The projected distance and velocity separation (the actual observed separation
maybe higher depending on the inclination angle) between the clouds is denoted
by $R$ and $\sigma_{v}$ respectively. From Torii et al. (2011), we adopt
$\sigma_{v}$ of 7.5 km s-1 and $R$ of 1–2 pc. From this, we estimate that the
mass, $M$ required to bind the clouds is $\sim$10${}^{4}\,M_{\odot}$, an order
of magnitude higher than both the stellar and cloud masses combined (Torii et
al., 2011). Hence the clouds can be considered to be gravitationally unbound.
Based on their velocity difference, they are expanding away from each other.
Note that Torii et al. (2011) arrived at a very similar conclusion, albeit for
a slightly higher distance estimate. Similarly, Torii et al. (2011) calculated
that the mechanical luminosity generated by the stellar wind of HD 164492A is
insufficient to cause to the observed cloud expansion. Therefore, we can
conclude similar to Torii et al. (2011) that the two molecular clouds are
expanding away from each other and are not gravitationally bound, and the
observed cloud velocities must be systemic. Based on the present direction of
the velocities and the cloud size, the two clouds must have collided $\sim$0.7
Myr ago (Torii et al., 2011), not accounting for projection effects. This
value may vary based on projection effects by a factor of 2.

The dynamical age ($t_{\textrm{dyn}}$) of the H II region encompassing the
Trifid nebula is estimated following;

$t_{\rm{dyn}}=\frac{4R_{s}}{7c_{s}}\left[\left(\frac{R_{{\rm H{\scriptsize
II}}}}{R_{s}}\right)^{(7/4)}-1\right].$ (4)

$R_{s}$ is the radius of the Strömgren sphere, $R_{\rm H{\scriptsize II}}$ is
the radius of the Trifid nebula, and $c_{s}$ is the sound speed. $R_{s}$ is
also calculated following eq. (1) in Tremblin et al. (2014). To calculate
$R_{s}$, a recombination coefficient of 3$\times 10^{-13}$ cm3 s-1 (adopting a
temperature of 8000 K), and a density of 3400 cm-3 is assumed from Tremblin et
al. (2014), as it best represents the solar neighbourhood. We also adopt the
ionising flux from Smith et al. (2002) for the central O7.5V star to be
5$\times 10^{48}$ photons s-1. Differences in the assumptions of mean density
can critically change the dynamical age (e.g. Tremblin et al. 2014), so we
calculate errors bars on our sample assuming a mean density range of
1000–10000 cm-3. From Fig. 1, we estimate $R_{\rm H{\scriptsize II}}\sim$2.5
pc, assuming a distance of 1250 pc. The resulting dynamical age is
$\sim$0.85${}^{+0.15}_{-0.45}$ Myr. This is larger than the estimate of 0.3
Myr arrived at by Cernicharo et al. (1998), assuming a distance of 1680 pc,
and density of 1000 cm-3. The value reported is in agreement with the
spectroscopic age of HD 164492A ($\sim$0.6 Myr).

####  5.1.2 Filamentary structure

Long-wavelength images of the Trifid Nebula taken from Herschel space
telescope display long filamentary structures. These structures can be
identified and extracted using various decomposition algorithms. To identify
the filamentary structures visible in the 500$\mu$m Herschel image of the
Trifid Nebula, we utilise the FilFinder algorithm from Koch & Rosolowsky
(2015). The ‘filaments’ are identified from a flattened and masked image to
remove bright compact sources, and pick out the filamentary structure. The
size threshold for masking chosen was experimented with to better highlight
the small scale and faint regions when visually comparing to the flattened
image, and a value of 600 was used. The final pruned skeletons representing
the ‘filaments’ are shown in Fig. 13. From the image, we identified two
filaments in the Trifid nebula. The two structures do not correspond visually
to the optical nebulosity, and the three visible dust lanes. The filaments run
along the visible emission in the mIR emission, but do not intersect.

A clear distinction can be made between the two filaments in terms of the
young stellar population. Along the northern filament lies the bulk of the
Class II YSO and PMS population detected, while only three fall along the
southern filament. The southern filament, especially towards the south west
edge falls across a clump of cold cores identified by Lefloch et al. (2002),
and a grouping of Class I YSOs. This suggests that stars forming along the
southern filament must be younger than the bulk of the stellar population
identified lying along the northern filament given the relative ratio of PMS
stars and Class II YSOs to Class I YSOs and cold cores.

####  5.1.3 Spatial distribution and star formation chronology of YSOs and
accreting PMS stars

The mean ages of the Class II YSOs are assumed to be similar or slightly
younger than the accreting PMS stars, whose median age is 1.5 Myr. The Class I
YSOs are considered to be much younger (Harvey et al., 2006), around $\sim$0.5
Myr. The dynamical age of the region ($\sim$0.85 Myr) is larger than the
estimated age of the Class II YSOs, with the former specifying a lower limit
to the age of the region. Further precise estimates of the reddening and
distance may verify the ages of the accreting PMS stars and Class II YSOs (see
Section 4.4). An additional timescale in consideration here is the crossing
time. Following Adamo et al. (2020), and assuming a total stellar mass of 500
$M_{\odot}$ based on Torii et al. (2011), the crossing time estimated is
$\sim$2.7 Myr. This is 2$\sigma$ larger than the median age of the accreting
PMS stars. Future spectroscopy of cluster members, that can estimate the
radial velocity spread of the region is essential to confirm the crossing time
(Adamo et al., 2020).

Fig. 13 shows the positions of the accreting PMS stars, YSOs, and cold cores
detected by Tapia et al. (2018) in the Trifid nebula. They are expected to
represent an age sequence of oldest to youngest stars currently forming in M
20. These positions are overlaid on the H$\alpha$ VPHAS+ image, the Spitzer
8$\mu$m image, and the Herschel 350$\mu$m respectively for the reader to
relatively gauge the position of the H$\alpha$ nebulosity, and the PAH and IR
emission. Also shown are the 2 km s-1 and cloud C contours taken from Torii et
al. (2011). The approximate location of the central cavity is also given, and
corresponds to the complementary distribution region in Torii et al. (2017)
between $18^{h}02^{m}15^{s}-19^{s}$, and
$-23^{\arcdeg}03^{\arcmin}-01^{\arcmin}$. The boundaries of the molecular
clouds can be visualised by the respective cloud contours given in Fig. 13.

From Fig. 13 the positions of the accreting PMS stars and YSOs are seen to
correlate to the positions of the molecular clouds. The PMS stars and Class II
YSOs are preferentially found along the boundaries of the molecular clouds,
with none found within the boundaries of the central cavity. Class II YSOs are
also located along the edges of the dust lanes. Further, the ages of the stars
show no clear demarcation between old or young stars, besides a group of older
stars lying outside the displayed boundaries of the molecular clouds to the
south east. The absence of PMS stars near the centre is not decisive evidence,
as the PMS sample is incomplete spatially. However, also no Class II YSOs are
detected towards the centre of the Nebula. This can be considered as evidence
as the Spitzer photometry is not affected towards the central region when
compared to the optical photometry.

In contrast, the Class I YSOs are found favouring the dust lanes near the
centre of the nebula. In addition, a concentration of Class I YSOs and cold
cores are found towards the South West, outside clearly visible nebulosity in
H$\alpha$, but near nebulous filaments visible in the 8$\mu$m image. This
suggests that a newer generation of stars is being formed near the periphery,
where the cold cores, and the Class I sources are located. A single flat
source is detected towards the central region, and few Class I sources, and
the cold core TC1 are found near HD 164992A.

From the ages and positions of the PMS stars, YSOs, and cold core samples
compared to the colliding molecular clouds in Fig. 13, there is evidence
suggestive of CCC induced star formation. The locations of the stars are along
the edges of the two clouds with similar ages. This age is smaller than the
crossing timescale $\sim$2.7 Myr. No stars are found near the centre of the
two colliding clouds, suggestive of a cavity devoid of on-going star
formation. The PMS stars, Class II YSOs have ages $\sim$1.5 Myr. Note that a
small grouping of Class I YSOs, and cold cores to the south west of the
nebula, falling on the southern filament are younger than most of the stellar
population in the region, and are forming in a subsequent outburst of star
formation.

In summary the gravitational energy of the clouds is not sufficient to bind
them, and the heating from HD 164492A seems insufficient to cause the observed
expansion. The velocity differences between the two clouds is systematic.
Travelling backwards, they must have collided $\sim$0.7 Myr. The collision of
the small cloud created a cavity on the larger cloud where the two clouds
collide and now overlap in the projected space (Torii et al., 2017). Star
formation was triggered likely immediately after, and molecular material
across the two clouds collapsed to form the currently visible young stellar
population roughly a million years ago. The positions and ages of YSOs and PMS
stars lend credence to the hypothesis that they may have formed after CCC. The
Class II YSOs and PMS stars ($\sim$1.5 Myr) are located preferentially towards
the edges of the clouds, and not found in the centre. Class I YSOs ($\lesssim
0.5$ Myr) are found more concentrated around the centre of the clouds. The
stars on either ends of the molecular cloud are around the same age, which is
less than the crossing time of $\sim$2.7 Myr.

A tightness in the dynamical age (0.85 Myr) and the age of the accreting PMS
stars may be alleviated if there are precise reddening, and distance
estimates. Further spectroscopic observations are needed to confirm this. Such
spectroscopic observations can also potentially identify discless/non-
accreting members. While the ages of such members are in keeping with the
general age of the region they reside in, they are often spatially less
concentrated (although occupying similar spatial distributions) compared to
their more active brethren.

##  6 Conclusions

  1. 1.

We have identified 46 Class II, 9 Flat, and 6 Class I YSOs in the Trifid
Nebula on the basis of their infrared SED slope, $\alpha$, and positions in
the mIR colour-colour plane. No Class III YSOs are included in our analysis.

  2. 2.

41 accreting PMS stars were identified on the basis of H$\alpha$ excess
against $ri$H$\alpha$ photometry. The identified stars have proper motions
indicating they likely belong to the Trifid Nebula. 90% of the accreting PMS
stars have nIR colours indicating the presence of a circumstellar dust disc.

  3. 3.

The distance to the Trifid Nebula measured from Gaia EDR3 parallaxes of a
subset of cluster members is around 1250 pc, placing the region in the
Sagittarius arm of the Milky Way.

  4. 4.

Based on the positions and ages of the young stellar population in the region,
it is likely that star formation was triggered by the collision of two clouds
$\sim$1 Myr ago, leading to the currently observed distribution.

V.M.K. acknowledges funding from CONICYT Programa de Astronomia Fondo Gemini-
Conicyt No 32RF180005. The work of V.M.K. is supported by NOIRLab, which is
managed by the Association of Universities for Research in Astronomy (AURA)
under a cooperative agreement with the National Science Foundation. V.M.K.
thanks K. Torii for kindly providing the molecular cloud cubes, and the
anonymous referee for detailed comments which helped improve this paper. This
work is based in part on observations made with the Spitzer Space Telescope,
which is operated by the Jet Propulsion Laboratory, California Institute of
Technology under a contract with NASA. This work is based in part on data
obtained as part of the UKIRT Infrared Deep Sky Survey. Based in part on
observations made with ESO Telescopes at the La Silla or Paranal Observatories
under programme ID(s) 177.D-3023(B), 177.D-3023(C), 177.D-3023(D),
177.D-3023(E). This work presents in part results from the European Space
Agency (ESA) space mission Gaia. Gaia data are being processed by the Gaia
Data Processing and Analysis Consortium (DPAC). Funding for the DPAC is
provided by national institutions, in particular the institutions
participating in the Gaia MultiLateral Agreement (MLA). Herschel is an ESA
space observatory with science instruments provided by European-led Principal
Investigator consortia and with important participation from NASA.

## References

  * Adamo et al. (2020) Adamo, A., Zeidler, P., Kruijssen, J. M. D., et al. 2020, Space Sci. Rev., 216, 69, doi: 10.1007/s11214-020-00690-x
  * Allen et al. (2004) Allen, L. E., Calvet, N., D’Alessio, P., et al. 2004, ApJS, 154, 363, doi: 10.1086/422715
  * Bailer-Jones et al. (2018) Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Mantelet, G., & Andrae, R. 2018, AJ, 156, 58, doi: 10.3847/1538-3881/aacb21
  * Barrado y Navascués & Martín (2003) Barrado y Navascués, D., & Martín, E. L. 2003, AJ, 126, 2997, doi: 10.1086/379673
  * Bressan et al. (2012) Bressan, A., Marigo, P., Girardi, L., et al. 2012, MNRAS, 427, 127, doi: 10.1111/j.1365-2966.2012.21948.x
  * Cambrésy et al. (2011) Cambrésy, L., Rho, J., Marshall, D. J., & Reach, W. T. 2011, A&A, 527, A141, doi: 10.1051/0004-6361/201015863
  * Cernicharo et al. (1998) Cernicharo, J., Lefloch, B., Cox, P., et al. 1998, Science, 282, 462, doi: 10.1126/science.282.5388.462
  * Churchwell et al. (2009) Churchwell, E., Babler, B. L., Meade, M. R., et al. 2009, PASP, 121, 213, doi: 10.1086/597811
  * Cutri et al. (2003) Cutri, R. M., Skrutskie, M. F., van Dyk, S., et al. 2003, 2MASS All Sky Catalog of point sources.
  * D’Alessio et al. (2001) D’Alessio, P., Calvet, N., & Hartmann, L. 2001, ApJ, 553, 321, doi: 10.1086/320655
  * De Marchi et al. (2010) De Marchi, G., Panagia, N., & Romaniello, M. 2010, ApJ, 715, 1, doi: 10.1088/0004-637X/715/1/1
  * Dotter et al. (2008) Dotter, A., Chaboyer, B., Jevremović, D., et al. 2008, ApJS, 178, 89, doi: 10.1086/589654
  * Drew et al. (2014) Drew, J. E., Gonzalez-Solares, E., Greimel, R., et al. 2014, MNRAS, 440, 2036, doi: 10.1093/mnras/stu394
  * Feigelson et al. (2013) Feigelson, E. D., Townsley, L. K., Broos, P. S., et al. 2013, ApJS, 209, 26, doi: 10.1088/0067-0049/209/2/26
  * Fukui et al. (2020) Fukui, Y., Habe, A., Inoue, T., Enokiya, R., & Tachihara, K. 2020, PASJ, doi: 10.1093/pasj/psaa103
  * Gaia Collaboration et al. (2020) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2020, arXiv e-prints, arXiv:2012.01533.  https://arxiv.org/abs/2012.01533
  * Gullbring et al. (1998) Gullbring, E., Hartmann, L., Briceno, C., & Calvet, N. 1998, ApJ, 492, 323, doi: 10.1086/305032
  * Gutermuth et al. (2009) Gutermuth, R. A., Megeath, S. T., Myers, P. C., et al. 2009, ApJS, 184, 18, doi: 10.1088/0067-0049/184/1/18
  * Harvey et al. (2006) Harvey, P. M., Chapman, N., Lai, S.-P., & Evans, Neal J., I. e. a. 2006, ApJ, 644, 307, doi: 10.1086/503520
  * Indebetouw et al. (2005) Indebetouw, R., Mathis, J. S., Babler, B. L., et al. 2005, ApJ, 619, 931, doi: 10.1086/426679
  * Kalari (2019) Kalari, V. M. 2019, MNRAS, 484, 5102, doi: 10.1093/mnras/stz250
  * Kalari et al. (2015) Kalari, V. M., Vink, J. S., Drew, J. E., et al. 2015, MNRAS, 453, 1026.  https://arxiv.org/abs/1507.06786
  * Koch & Rosolowsky (2015) Koch, E. W., & Rosolowsky, E. W. 2015, MNRAS, 452, 3435, doi: 10.1093/mnras/stv1521
  * Lefloch et al. (2002) Lefloch, B., Cernicharo, J., Rodríguez, L. F., et al. 2002, ApJ, 581, 335, doi: 10.1086/344049
  * Loren (1976) Loren, R. B. 1976, ApJ, 209, 466, doi: 10.1086/154741
  * Lucas et al. (2008) Lucas, P. W., Hoare, M. G., Longmore, A., et al. 2008, MNRAS, 391, 136, doi: 10.1111/j.1365-2966.2008.13924.x
  * McKee & Ostriker (2007) McKee, C. F., & Ostriker, E. C. 2007, ARA&A, 45, 565, doi: 10.1146/annurev.astro.45.051806.110602
  * Megeath et al. (2004) Megeath, S. T., Allen, L. E., Gutermuth, R. A., et al. 2004, ApJS, 154, 367, doi: 10.1086/422823
  * Meyer et al. (1997) Meyer, M. R., Calvet, N., & Hillenbrand, L. A. 1997, AJ, 114, 288, doi: 10.1086/118474
  * Pecaut & Mamajek (2013) Pecaut, M. J., & Mamajek, E. E. 2013, ApJS, 208, 9, doi: 10.1088/0067-0049/208/1/9
  * Petit et al. (2019) Petit, V., Wade, G. A., Schneider, F. R. N., et al. 2019, MNRAS, 489, 5669, doi: 10.1093/mnras/stz2469
  * Pickles (1998) Pickles, A. J. 1998, PASP, 110, 863, doi: 10.1086/316197
  * Rho et al. (2001) Rho, J., Corcoran, M. F., Chu, Y.-H., & Reach, W. T. 2001, ApJ, 562, 446, doi: 10.1086/323053
  * Rho et al. (2008) Rho, J., Lefloch, B., Reach, W. T., & Cernicharo, J. 2008, M20: Star Formation in a Young HII Region, ed. B. Reipurth, Vol. 5, 509
  * Rho et al. (2006) Rho, J., Reach, W. T., Lefloch, B., & Fazio, G. G. 2006, ApJ, 643, 965, doi: 10.1086/503245
  * Robitaille & Bressert (2012) Robitaille, T., & Bressert, E. 2012, APLpy: Astronomical Plotting Library in Python.  http://ascl.net/1208.017
  * Robitaille et al. (2006) Robitaille, T. P., Whitney, B. A., Indebetouw, R., Wood, K., & Denzmore, P. 2006, ApJS, 167, 256, doi: 10.1086/508424
  * Siess et al. (2000) Siess, L., Dufour, E., & Forestini, M. 2000, A&A, 358, 593
  * Smith et al. (2002) Smith, L. J., Norris, R. P. F., & Crowther, P. A. 2002, MNRAS, 337, 1309, doi: 10.1046/j.1365-8711.2002.06042.x
  * Tapia et al. (2018) Tapia, M., Persi, P., Román-Zúñiga, C., et al. 2018, MNRAS, 475, 3029, doi: 10.1093/mnras/sty048
  * Torii et al. (2011) Torii, K., Enokiya, R., Sano, H., et al. 2011, ApJ, 738, 46, doi: 10.1088/0004-637X/738/1/46
  * Torii et al. (2017) Torii, K., Hattori, Y., Hasegawa, K., et al. 2017, ApJ, 835, 142, doi: 10.3847/1538-4357/835/2/142
  * Tremblin et al. (2014) Tremblin, P., Anderson, L. D., Didelon, P., et al. 2014, A&A, 568, A4, doi: 10.1051/0004-6361/201423959
  * Urquhart et al. (2014) Urquhart, J. S., Csengeri, T., Wyrowski, F., et al. 2014, A&A, 568, A41, doi: 10.1051/0004-6361/201424126
  * Wade et al. (2017) Wade, G. A., Shultz, M., Sikora, J., et al. 2017, MNRAS, 465, 2517, doi: 10.1093/mnras/stw2799
  * Yusef-Zadeh et al. (2005) Yusef-Zadeh, F., Biretta, J., & Geballe, T. R. 2005, AJ, 130, 1171, doi: 10.1086/432095
  * Yusef-Zadeh et al. (2000) Yusef-Zadeh, F., Shure, M., Wardle, M., & Kassim, N. 2000, ApJ, 540, 842, doi: 10.1086/309352