A potential new phase of massive star formation

L. Bonne; N. Peretto; A. Duarte-Cabral; A. Schmiedeke; N. Schneider; S. Bontemps; A. Whitworth

doi:10.1051/0004-6361/202142154

Home

All issues

Volume 665 (September 2022)

A&A, 665 (2022) A22

Full HTML

Open Access

Issue		A&A Volume 665, September 2022


Article Number		A22
Number of page(s)		22
Section		Interstellar and circumstellar matter
DOI		https://doi.org/10.1051/0004-6361/202142154
Published online		06 September 2022

A&A 665, A22 (2022)

A luminous outflow cavity centred on an infrared quiet core

L. Bonne¹^,2^,3, N. Peretto¹, A. Duarte-Cabral¹, A. Schmiedeke⁴, N. Schneider⁵, S. Bontemps² and A. Whitworth¹

¹ School of Physics & Astronomy, Cardiff University, Queens Buildings, The parade, Cardiff CF24 3AA UK
e-mail: lbonne@usra.edu
² Laboratoire d’Astrophysique de Bordeaux, Université de Bordeaux, CNRS, B18N, allée Geoffrey Saint-Hilaire, 33615 Pessac, France
³ SOFIA Science Center, USRA, NASA Ames Research Center, Moffett Field, CA 94045 USA
⁴ Max-Planck-Institut für extraterrestrische Physik, Gießenbachstraße 1, 85748 Garching, Germany
⁵ I. Physik. Institut, University of Cologne, Zülpicher Str. 77, 50937 Cologne, Germany

Received: 3 September 2021
Accepted: 3 May 2022

Abstract

Context. Due to the sparsity and rapid evolution of high-mass stars, a detailed picture of the evolutionary sequence of massive protostellar objects still remains to be drawn. Some of the early phases of their formation are so short that only a select number of objects throughout the Milky Way currently find themselves spending time in those phases.

Aims. Star-forming regions going through the shortest stages of massive star formation present different observational characteristics than most regions. By studying the dust continuum and line emission of such unusual clouds, one might be able to set strong constraints on the evolution of massive protostellar objects.

Methods. We present a detailed analysis of the G345.88-1.10 hub filament system, which is a newly discovered star-forming cloud that hosts an unusually bright bipolar infrared nebulosity at its centre. We used archival continuum observations from Berschel, WISE, Spitzer, 2MASS, and SUMSS in order to fully characterise the morphology and spectral energy distribution of the region. We further made use of APEX ¹²CO(2–1), ¹³CO(2–1), C¹⁸O(2–1), and H30α observations to investigate the presence of outflows and map the kinematics of the cloud. Finally, we performed RADMC-3D radiative transfer calculations to constrain the physical origin of the central nebulosity.

Results. At a distance of 2.26_-0.21^+0.30 kpc, G345.88-1.10 exhibits a network of parsec-long converging filaments. At the junction of these filaments lie four infrared-quiet fragments. The fragment H1 is the densest one (with M = 210 M_⊙, R_eff = 0.14 pc) and sits right at the centre of a wide (opening angle of ~90 ± 15°) bipolar nebulosity where the column density reaches local minima. The ¹²CO(2–1) observations of the region show that these infrared-bright cavities are spatially associated with a powerful molecular outflow that is centred on the H1 fragment. Negligible radio continuum and no H30α emission is detected towards the cavities, seemingly excluding the idea that ionising radiation drives the evolution of the cavities. Radiative transfer calculations of an embedded source surrounded by a disc and/or a dense core are unable to reproduce the observed combination of a low-luminosity (≲500 L_⊙) central source and a surrounding high-luminosity (~4000 L_⊙) mid-infrared-bright bipolar cavity. This suggests that radiative heating from a central protostar cannot be responsible for the illumination of the outflow cavities.

Conclusions. This is, to our knowledge, the first reported object of this type. The rarity of objects such as G345.88-1.10 is likely related to a very short phase in the massive star and/or cluster formation process that has been unidentified thus far. We discuss whether mechanical energy deposition by one episode or successive episodes of powerful mass accretion in a collapsing hub might explain the observations. While promising in some aspects, a fully coherent scenario that explains the presence of a luminous bipolar cavity centred on an infrared-dark fragment remains elusive at this point.

Key words: stars: protostars / ISM: jets and outflows / stars: formation / ISM: kinematics and dynamics / ISM: clouds

© L. Bonne et al. 2022

Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe-to-Open model. Subscribe to A&A to support open access publication.

1 Introduction

The series of processes that lead to the transport of ambient molecular gas down to the surface of young accreting protostars is not fully understood yet (Li et al. 2014). Questions related to the nature of the initial instability, the star formation timescale and the morphology and size of the mass reservoir are examples of areas in which progress needs to be made. The study of nearby (d < 500 pc) star-forming regions presents a picture in which prestellar cores, that is the progenitors of individual stars or small stellar systems (Ward-Thompson et al. 2007; di Francesco et al. 2007), form as the result of the fragmentation of gravitationally unstable filaments (André et al. 2014). These cores then collapse to form a protostar at their centres. Conservation of angular momentum means that a circumstellar disc is formed alongside the protostar. The size and time at which such discs first appear, along with the exact role they play in the accretion process, are still a matter of active debate (e.g. Li et al. 2014; Hartmann et al. 2016). As far as low-mass star formation is concerned, most of the stellar mass is thought to be accreted during the first (i.e. Class 0) stages (André et al. 2000). It is important to note that there is now plenty of evidence that accretion is not a steady process but rather episodic (e.g. Hartmann et al. 2016). A consequence of the mass accretion activity is the formation of high-velocity collimated jets that evacuate a part of the angular momentum from the core during protostellar accretion. This jet entrains with it ambient molecular gas that is seen as a molecular protostellar outflow (Bachiller 1996; Bally 2016). The presence of a protostellar outflow thus provides a probe of protostellar mass accretion. These outflows might also contribute to the self-regulation of star formation within molecular clouds by injecting turbulence and slowing down star formation activity (e.g. Arce et al. 2011; Duarte-Cabral et al. 2012).

The results presented above have mostly been obtained by studying low-mass star-forming regions. It is however debated whether this applies to the formation of the rarer and faster evolving high-mass stars (Zinnecker & Yorke 2007; Motte et al. 2018). For instance, despite extensive searches in the Galactic plane, no significant population of massive (M > 100 M_⊙, R ≤ 0.05 pc) prestellar cores has been found (Motte et al. 2007, 2010; Russeil et al. 2010; Louvet et al. 2019; Sanhueza et al. 2019). Instead, a growing body of evidence is showing that the mass of massive protostellar sources increases, at least initially, with time (Tigé et al. 2017; Peretto et al. 2020; Rigby et al. 2021), most likely as the result of the global collapse of their parent clumps (Peretto et al. 2006, 2013, 2014; Schneider et al. 2010, 2015a; Kirk et al. 2013; Csengeri et al. 2017a,b; Watkins et al. 2019; Jackson et al. 2019). The presence of hub filament systems, i.e. networks of converging filaments, might be a tracer of such collapse. In fact, it has been shown that the centre of such systems represent ideal locations for the formation of massive cores (e.g. Schneider et al. 2012; Williams et al. 2018), with hubs being more efficient at concentrating a large fraction of their mass within the most massive cores (Anderson et al. 2021). The large-scale kinematics of the parent cloud around massive cores is therefore key to better understand high-mass star formation. A direct link between the mass infall rate of the prototypical hub SDC335 and its mass outflow rate even suggests a continuity of mass transfer from parsec-scales down to core scale (Avison et al. 2021).

Another major difference between the formation of low-mass and high-mass stars is their relative Kelvin-Helmoltz time scale, i.e. the time it takes for pre-main sequence stars to radiate away its thermal energy. For low-mass stars, this time is larger than the lifetime of pre-main sequence stars, allowing low-mass stars to start the ignition of hydrogen burning after it finished accreting. This is not the case for high-mass stars (>10 M_⊙) which have shorter Kelvin-Helmoltz timescales than the accretion timescale. As a result, these stars start burning hydrogen while still deeply embedded within their parent molecular cloud. This used to be a major hurdle for our understanding of high-mass star formation since the large radiation pressure from massive stars should halt accretion (Larson & Starrfield 1971; Wolfire & Cassinelli 1987). However, this is resolved when considering non-spherically symmetric accretion flows, such as through an accretion disc. Accretion discs provide a low-density channel for protostellar feedback to escape, the so-called flashlight effect, preventing the reversal of the mass inflow (e.g. Yorke & Bodenheimer 1999; Yorke & Sonnhalter 2002; Kuiper et al. 2010). With the advent of (sub)millimetre high-angular resolution observatories such as ALMA and NOEMA, an increasing number of massive discs/torus have been found (Johnston et al. 2015; Beltrán & de Wit 2016; Ilee et al. 2016; Cesaroni et al. 2017; Csengeri et al. 2018). However, these discs are often fragmented, a possible consequence of the high density of these structures (Beuther et al. 2018; Ilee et al. 2018; Ahmadi et al. 2018). Like in the low-mass case, it has been proposed that high-mass protostars undergo episodic accretion bursts (e.g. Carrasco-González et al. 2015; Stecklum et al. 2021; Hunter et al. 2021).

In order to make further progress in our understanding of the earliest stages of massive star formation, two approaches can be considered. The first one consists in studying a large sample of sources in order to get statistics on a set of properties. The second one consists in performing a detailed analysis of an individual source that, for one reason or another, is of particular interest. It is this second approach that we are taking in this paper. We present here the first analysis of the G345.88-1.10 star-forming cloud (see Fig. 1), a hub filament system with a bright infrared bipolar nebulosity at the centre. G345.88-1.10 has been serendipitously discovered, and, as of today, no literature exists on it. As discussed in this paper, we believe that the study of G345.88-1.10 will shed new light on the mechanisms linked to the very first phases of high-mass star formation. In Sect. 2 we describe the observational data we used in the present study. Section 3 presents the hub characteristics. Section 4 describes the morphology and spectral energy distribution (SED) of the bipolar nebulosity. In Sect. 5 we identify and analyse the properties of the central molecular outflow, while in Sect. 6 we investigate the possibility whether the bipolar nebulosity could be the result of ionising radiation. In Sect. 7 the radiative transfer (RT) calculations are presented and compared with the observed Herschel continuum emission. In Sect. 8, we bring all the results together, and discuss their implications for massive star/cluster formation. Finally, Sect. 9 presents a summary of the paper and future prospects.

2 Observational Data

As a result of the relatively large galactic latitude of G345.88-1.10 (|b| > 1) a few key Galactic plane surveys performed in the past 20 yr, such as MIPSGAL (Carey et al. 2009) or more recently SEDIGISM (Schuller et al. 2017, 2021), have not observed it. This could explain why it did not attract attention yet. In this paper, we study G345.88-1.10 across a number of wavelengths, using dust continuum data from Herschel, Spitzer, WISE, and new spectral line observations obtained with the Atacama Pathfinder Experiment telescope (APEX).

2.1 Herschel

The Herschel data of G345.88-1.10 is provided by the Hi-GAL survey which is dedicated to the far-infrared dust continuum imaging of the Galactic plane, and was obtained from the HiGAL archive¹ (Molinari et al. 2010, 2016; Traficante et al. 2011). The photometric maps at 70 and 160 μm, obtained with the PACS instrument (Poglitsch et al. 2010), have an angular resolution of 8.4″ and 13.5″, respectively. The maps at 250, 350 and 500 μm, obtained with the SPIRE instrument (Griffin et al. 2010), have a resolution of 18.2″, 24.9″ and 36.3″, respectively.

2.2 WISE, Spitzer and 2MASS IR surveys

The WISE all-sky survey provides data at 3.4, 4.6, 12 and 22 μm for G345.88-1.10 with an angular resolution of 6.1″, 6.4″, 6.5″ and 12.0″, respectively (Wright et al. 2010).

Although G345.88-1.10 has not been observed with the MIPS instrument on the Spitzer space telescope, the source has been covered by the IRAC instrument within the GLIMPSE survey (Benjamin et al. 2003), providing images at 3.6, 4.5, 5.8 and 8 μm at an angular resolution of 1.66, 1.72, 1.88 and 1.98″, respectively (Fazio et al. 2004).

Lastly, the 2MASS survey provides a map in the J, H and K-band for G345.88-1.10 with an angular resolution of ~2″ (Skrutskie et al. 2006).

Fig. 1

Left: H₂ column density image of G345.88-1.10 on a log scale at a spatial resolution of 18″ derived from Herschel data at 160 and 250 μm. The hub morphology has been highlighted by the skeleton of the filament network present in the region. Right: corresponding Herschel dust temperature map. The black box indicates the region of interest for this paper which is presented in more detail in Fig. 3.

2.3 APEX

In march 2016, the ¹²CO(2–1), ¹³CO(2–1) and C¹⁸O(2–1) lines were mapped around G345.88-1.10 using the APEX-1 receiver which is installed on the APEX telescope (Güsten et al. 2006). The frequencies (and wavelengths) of these rotational transitions are 230.54 (1.301), 220.40 (1.361) and 219.56 GHz (1.365 mm), respectively. The setup that covers ¹²CO(2–1) simultaneously also covers the H30α hydrogen recombination line at 231.90 GHz, which is a tracer of Hii regions. This data has an angular resolution of ~28″. The ¹³CO(2–1) and C¹⁸O(2–1) lines were mapped simultaneously by tuning the receiver to an intermediate frequency. All lines were observed with a precipitable water vapor (pwv) between 0.8 and 1.3 mm and the data was converted to the main beam temperature (T_mb) using a main beam efficiency η_mb = 0.81 (Vassilev et al. 2008). The final maps were produced using the GILDAS CLASS software² with a channel width of 0.3 km s⁻¹ for all lines. The noise, estimated as the line standard deviation within signal-free channels, is ~0.2 K for ¹²CO(2–1) and H30α. For ¹³CO(2–1) and C¹⁸O(2−1), the typical noise over the map is ~0.15 K. The observations were carried out using the on-the-fly mode centred on α₍₂₀₀₀₎ = 17^h 11^m43.9^s, δ₍₂₀₀₀₎ = −41° 18′25.5″. The ¹²CO(2–1) and H30α map covers an area on the sky of 4′ × 4′ while the ¹³CO(2–1) and C¹⁸O(2–1) maps cover a smaller area of 2′ × 2′.

3 G345.88-1.10: A Nearby Hub Filament System

3.1 The Distance of G345.88-1.10

As this is the first paper on G345.88-1.10, no published distance exists. To estimate it, we use the galactic rotation model from Reid et al. (2014). 0ur APEX data show (see Fig. 2) that the velocity of G345.88-1.10 is centred at −21 km s⁻¹. This results in a heliocentric kinematic near distance of ${2.26}_{- 0.21}^{+ 0.30}$ $2.26_{ - 0.21}^{ + 0.30}$ kpc, where the uncertainties are given by the Monte Carlo calculations³ presented in Wenger et al. (2018). This distance estimate places G345.88-1.10 in the Scrutum galactic arm (e.g. Pineda et al. 2013), at a galactocentric distance of 6.1 kpc. The far distance gives ${13.96}_{- 0.48}^{+ 0.24}$ $13.96_{ - 0.48}^{ + 0.24}$ pc. However, this far distance seems unlikely to us for a couple of reasons. First, it would imply that the bolometric luminosity of the cavity exceeds 10⁵ L_⊙ (see Sect. 4.2 for the L_bol estimate). Galactic sources that bright are generally associated with bright radio continuum emission (>100 mJy) which we do not observe in G345.88-1.10 (see Sect. 6). We note that a couple of luminous (L_bol > 10⁵ L_⊙) sources were found with a radio continuum flux <100 mJy (Lumsden et al. 2013), but only a small fraction of luminous sources satisfies this condition. Second, if placed at the far distance, then the two clouds at −10 km s⁻¹ and −6 km s⁻¹ that we identify in Fig. 2 would be foreground clouds. These two clouds overlap in velocity with the observed C0 high-velocity wings in Fig. 2. Given that outflows in dense molecular regions typically have higher excitation temperatures than the bulk of diffuse molecular gas, at least a part of these clouds would likely be seen in absorption rather than emission. This argument is particularly valid for G345.88-1.10 as we show in Sect. 5 that the outflow of G345.88-1.10 is located in the warm nebulosities shown in Fig. 1. Lastly, we note that assuming the far distance for the region implies that the region is located 270 pc above the Galactic plane. This is at the edge of the thin disk in the Milky Way. As the region is forming high-mass stars, see later on, this is a fairly high distance from the Galactic plane. For the remainder of this paper, we therefore select the close distance solution for G345.88-1.10.

Fig. 2

¹²CO(2–1), ¹³CO(2–1), and C¹⁸O(2–1) emission averaged over the extent of the hub.

Fig. 3

Left: Berschel 70 μm image of the central region/hub of G345.88-1.10, showing the presence of a bipolar nebulosity N1 and N2. The H₂ column density contours are overplotted in white starting at $N_{H_{2}} = 2 \times 10^{22}$ ${N_{{{\rm{H}}_2}}} = 2 \times {10^{22}}$ cm⁻² with increments of 2 × 10²² cm⁻². The 4 fragments that were extracted from the column density map with the dendrogram technique are indicated by the white crosses together with the two nebulosities. Right: image of the Berschel dust temperature centred on the infrared nebulosity. All symbols and notations are identical to those presented in the left-hand-side panel.

3.2 H₂ Column Density and Dust Temperature Mapping

Figure 1 presents the H₂ column density and dust temperature images of G345.88-1.10 at 18″ resolution. These have been obtained using the method described in Peretto et al. (2016). This method only uses the 160 and 250 μm data from Berschel. First, we convolved the 160 μm image to the resolution of the 250 μm image using a Gaussian kernel of $F W H M = \sqrt{θ_{250}^{2} - θ_{160}^{2}} = 13. 4^{″}$ $FWHM = \sqrt {\theta _{250}^2 - \theta _{160}^2} = 13.4''$ . Then, the temperature at every pixel was determined from the ratio of the 160 and 250 μm flux density using $R_{160 / 250} = \frac{S_{160}}{S_{250}} = \frac{B_{v 160} (T_{d})}{B_{v 250} (T_{d})} {(\frac{250}{160})}^{β}$ ${R_{{{160} \mathord{\left/ {\vphantom {{160} {250}}} \right. \kern-\nulldelimiterspace} {250}}}} = {{{S_{160}}} \over {{S_{250}}}} = {{{B_{v160}}\left( {{T_d}} \right)} \over {{B_{v250}}\left( {{T_d}} \right)}}{\left( {{{250} \over {160}}} \right)^\beta }$ (1)

where S_λ is the flux density, B_ν the Planck function, β the index of the dust opacity law and T_d the dust temperature. For β a value of 2 was used (Sadavoy et al. 2013). Since Eq. (1) has no analytical solution, Brent’s method⁴ was used to find a zero point within the 10–100 K temperature range. For all pixels, a solution was found in the 13–30 K range. From the resulting dust temperature map, the H₂ column density image can then be obtained via the following equation: $N_{H_{2}} = S_{v_{250}} / [B_{v_{250}} {(T_{d})}_{κ_{v_{250}}} {μm}_{H}]$ ${N_{{{\rm{H}}_2}}} = {{{S_{{v_{250}}}}} \mathord{\left/ {\vphantom {{{S_{{v_{250}}}}} {[{B_{{v_{250}}}}{{\left( {{T_d}} \right)}_{{\kappa _{{v_{250}}}}}}{\rm{\mu }}{{\rm{m}}_H}]}}} \right. \kern-\nulldelimiterspace} {[{B_{{v_{250}}}}{{\left( {{T_d}} \right)}_{{\kappa _{{v_{250}}}}}}{\rm{\mu }}{{\rm{m}}_H}]}}$ (2)

where $κ_{v_{250}} = 0.012$ ${\kappa _{{v_{250}}}} = 0.012$ cm² g⁻¹ is the specific dust opacity at 250 μm that already incorporates a dust to gas mass ratio of 1% (Hildebrand 1983), μ = 2.33 is the average molecular weight and m_H is the mass of atomic hydrogen.

As a safety-check, we also produced the temperature and column density maps by performing a pixel-by-pixel SED fit to the 160, 250, 350 and 500 μm maps, at an angular resolution of 36″ (e.g. Peretto et al. 2010; Battersby et al. 2011; Hill et al. 2011). The resulting maps were very similar to those presented in Fig. 1, but at a lower angular resolution. Furthermore, the calculation of the total gas mass within the central parsec of G345.88-1.10 using both column density maps provided similar values within 10% of each other, i.e. 2900 M_⊙ for the ratio method and 2600 M_⊙ for the SED method. In the remainder of this paper, we consider the column density and temperature map produced by the ratio method of Peretto et al. (2016).

3.3 Identifying Filaments and Fragments

Figure 1 shows two very interesting features. First, there is a clear hub filament system with a network of a few parsec long filaments converging towards a high column density region. In order to highlight this structure, we have performed a skeleton extraction using the second derivative method presented in, for example, Schisano et al. (2014), Williams et al. (2018), Watkins et al. (2019), and Orkisz et al. (2019). The second striking feature is the presence of two warm regions right at the centre of G345.88-1.10, on either side of the highest column density peak.

In order to identify sub-structures within the column density map we used the astrodendro package⁵ (Rosolowsky et al. 2008). Astrodendro segments the input image by constructing a dendrogram tree. The most compact non-fragmented structures of this tree are referred to as leaves in the dendrogram terminology, and fragments in the remainder of this paper. To be identified as a fragment that is spatially resolved, the leaf has to cover at least two times the beamsize at 250 μm. This limits the resolution of the dendrogram to an effective radius of ~0.1 pc. The tree starts at $N_{H_{2}} = 10^{22}$ ${N_{{{\rm{H}}_2}}} = {10^{22}}$ cm⁻², to focus on the dense gas only, and is constructed using steps of $N_{H_{2}} = 2.5 \times 10^{21}$ ${N_{{{\rm{H}}_2}}} = 2.5 \times {10^{21}}$ cm⁻². With these parameters, astrodendro identifies 4 fragments in the hub. These are labelled as H1, H2, H3 and H4 in Fig. 3.

The fragment with the highest column density (H1) has a mass of $210_{- 30}^{+ 50}$ $210_{ - 30}^{ + 50}$ M_⊙ within a radius of 0.14 pc. The coordinates, mass, effective radius and average column density of the identified fragments are given in Table 1. The other identified fragments (H2–H4) also have masses in the 100–300 M_⊙ range. While we are expecting these sources to be sub-fragmented on smaller scales, their masses suggest that G345.88-1.10 has the potential to form a rich stellar cluster.

Table 1

Identified fragment properties.

4 A Luminous Bipolar Nebulosity Embedded within G345.88-1.10

4.1 Morphology of the Mid-Infrared Emission

As mentioned in the introduction, the striking feature of G345.88-1.10 is its central mid-infrared bipolar nebulosity (see Fig. 3). The nebulosity is made of two lobes, N1 and N2, that are located on either side of the H1 fragment. Figure 3 also shows that, in addition to be warmer, the dust within the nebulosity reaches local H₂ column density minima. Some 70 μm emission was recently found in the bipolar cavities of S106 (Adams et al. 2015) and RCW 36 (Minier et al. 2013) which both contain one or two young O stars at their centre (Ellerbroek et al. 2013; Comerόn et al. 2018; Schneider et al. 2018). However, what is striking about the 70 μm emission towards G345.88-1.10 is the dark H1 fragment in which no clear 70 μm source is observed. If a (proto)star, or a group of (proto)stars, is responsible for powering the mid-infrared nebulosities then it must be deeply embedded within H1. This 70 μm darkness of the H1 fragment is in stark contrast to S106 and RCW 36, where the 70 μm emission at the location of the embedded source is an order of magnitude brighter than the cavities (Minier et al. 2013; Adams et al. 2015).

In order to check whether or not the 70 μm nebulosities might be powered by two independent young clusters, we studied the literature and inspected images at shorter wavelengths down to the J-band. After eye inspection of the VISTA images, Borissova et al. (2011) proposed that G345.88-1.10 hosts two small clusters that are spatially coincident with the mid-infrared nebulosities. For each cluster, they identified 10 potential members which is exactly the lower limit they set in their study to identify a cluster candidate. However, inspecting the 2MASS RGB image of the J-band, B-band and K-band in Fig. 4 we find no convincing indication of a stellar over-density compared to the field stars. In fact, in the systematic 2MASS search for star clusters by Froebrich et al. (2007) no cluster candidates were detected within the G345.88-1.10 region. Furthermore, to our knowledge, no cluster candidates were found towards G345.88-1.10 in the Gaia database either. At longer wavelengths (see Fig. 5) we do see some compact sources at 3.6 μm within the boundary of the nebulosities, however, once again, the compact source density is not larger than that of the background. Another argument against the presence of two star clusters is the resemblance of the two nebulosities, both in terms of size and brightness. This is visible by eye (see Fig. 5), but will be quantified further in Sect. 4.2. Note that it would also be a remarkable coincidence that both clusters would appear to form a bipolar structure around the fragment with the highest density within the hub. Overall, it seems unlikely that the two infrared nebulosities observed towards G345.88-1.10 are powered by stellar clusters.

Finally, looking at Fig. 3 in more detail, two relatively weak compact 70 μm sources can be spotted, one nearly coincident with the dust column density peak of H3, and the other one just at the boundary of H4. These two sources have emission counterparts at all wavelengths presented in Fig. 5, demonstrating that these are protostellar sources. In fact, of the four fragments identified, only the H2 fragment shows absolutely no association with mid-infrared emission.

Fig. 4

RGB image of the 2MASS J-band (in blue), H-band (in green) and K-band (in red). The white contours indicate the 70 μm continuum emission at 2000 and 8000 MJy sr⁻¹, to highlight the location of the full 70 μm nebulosities and the brightest regions of the 70 μm nebulosities, respectively. The 2MASS data shows no evident stellar overdensity in the cavities which would point to the presence of clusters.

Fig. 5

Mid-infrared images from IRAC (3.6, 4.5, 5.8 and 8 μm), WISE (22 μm) and Herschel (160, 250, 350, 500 μm) of the G345.88-1.10 hub with the overlaid contour indicating the area of the 70 μm nebulosities (here defined by being brighter than 2000 MJy sr⁻¹).

4.2 The Spectral Energy Distribution and Luminosity

In order to determine the respective luminosities of the nebulosities and the H1 fragment we construct the SEDs of N1, N2, and H1 separately. For that purpose, fluxes were obtained from the Herschel, WISE and IRAC photometric maps using aperture photometry (see Table 2, Appendix A). For the H1 fragment, no fluxes were determined at wavelengths shorter than 160 μm since it is not possible to disentangle the fragment from the nebulosities at these shorter wavelengths. For the nebulosities no fluxes are calculated at wavelengths above 250 μm as they are probably associated with the ambient cloud. The corresponding SEDs are presented in Fig. 6. The first feature one notices is the strong similarity between the shape and magnitude of the N1 and N2 SEDs. 0nly at wavelengths covered by Herschel, the southern nebulosity (N2) becomes slightly brighter. The luminosity of each source is determined by integrating over the entire wavelength range, from 3.6 to 160 μm for N1 and N2, and from 160 to 500 μm for H1. For the two nebulosities, we obtain luminosities of ${2.1}_{- 1.1}^{+ 1.2} \times 10^{3}$ $2.1_{ - 1.1}^{ + 1.2} \times {10^3}$ L_⊙ and ${2.3}_{- 1.2}^{+ 1.3} \times 10^{3}$ $2.3_{ - 1.2}^{ + 1.3} \times {10^3}$ L_⊙, respectively. The luminosity of the H1 fragment was found to be $51_{- 27}^{+ 29}$ $51_{ - 27}^{ + 29}$ L_⊙ only. The luminosity of the entire region (N1 + N2 + H1) is 4.4 × 10³ L_⊙, and thus dominated by the nebulosities. When considering the upper limit of 5 Jy at 70 μm for the H1 fragment, an upper limit of 100 L_⊙ is obtained for the H1 fragment, which is still a fraction (i.e. ~2%) of the total luminosity. For a standard protostellar core, the sum of the core and cavity luminosities is generally considered to be representative of the protostellar luminosity, but we show later that most of the flux from the nebulosities cannot be related to photon emission originating from an embedded object in H1.

5 Associated Protostellar Outflows

5.1 The Location of High-Velocity Wings

Being centred on the H1 fragment, the bipolar nebulosity at the centre of G345.88-1.10 is reminiscent of outflow cavities. The presence of outflows can be verified by the inspection of the ¹²CO(2–1) data. In fact, in Fig. 2, one can already notice the presence of high-velocity wings in the average ¹²CO(2–1) spectrum of the hub, a clear signature of outflowing gas. By integrating these ¹²CO(2–1) high velocity wings, the location of molecular outflows within G345.88-1.10 can be established. While one could use fixed velocity intervals for the integration of the ¹²CO(2–1) wing emission, the presence of velocity gradients within the bulk of the cloud gas would typically lead to large errors. In order to minimise such a source of error, the C¹⁸ O(2–1) centroid velocity $(υ_{c^{18}}_{_{O}})$ $\left( {{\upsilon _{{{\rm{c}}^{18}}}}_{_{\rm{O}}}} \right)$ and velocity dispersion $(σ_{c^{18}}_{_{O}})$ $\left( {{\sigma _{{{\rm{c}}^{18}}}}_{_{\rm{O}}}} \right)$ were first derived for each pixel (see Appendix C). After visual inspection, we defined the ¹²CO outflow velocity ranges as $υ_{c^{18}}_{_{O}} + 3.5 σ_{c^{18}}_{_{O}} \leq υ_{red} \leq 0$ ${\upsilon _{{{\rm{c}}^{18}}}}_{_{\rm{O}}} + 3.5{\sigma _{{{\rm{c}}^{18}}}}_{_{\rm{O}}} \le {\upsilon _{{\rm{red}}}} \le 0$ and $- 40 \leq υ_{blue} \leq υ_{c^{18}}_{_{O}} - 3.5 σ_{c^{18}}_{_{O}}$ $- 40\, \le {\upsilon _{{\rm{blue}}}} \le \,{\upsilon _{{{\rm{c}}^{18}}}}_{_{\rm{O}}} - 3.5{\sigma _{{{\rm{c}}^{18}}}}_{_{\rm{O}}}$ . Finally, as can be seen in Fig. 2, contaminating emission is present around v_lsr = −5 km s⁻¹ km s⁻¹ and v_lsr = −10 km s⁻¹. To avoid contribution from these background clouds, the emission between −6.5 km s⁻¹ to −4 km s⁻¹ and −11 km s⁻¹ to −9 km s⁻¹ is left out of the integration of the wings. The resulting integrated intensity maps are represented as contours in Fig. 7 (top) and show the presence of two unresolved bipolar outflows. The kinematics of these outflows are presented in more detail in Appendix B, showing no clear signs of further substructure. The brightest of the two bipolar outflows is centred on the H1 fragment, along an axis that is oriented in the same direction as the bipolar nebulosities. The outflow is the brightest towards the region where the nebulosities’ edges exhibit the largest H₂ column density. Since molecular outflows are considered to be ambient gas that is being entrained by the release of momentum along the jet axis (Gueth et al. 1998; Frank et al. 2014), this slight misalignment is to be expected. This is also in agreement with high angular resolution observations of protostellar molecular outflows which show that the low-velocity ¹²CO outflow is best detected towards the edges of the outflow cavity (e.g. Tafalla 2017; Tabone et al. 2017; Louvet et al. 2018; de Valon et al. 2020). Altogether, the association between the identified outflow and the bipolar nebulosity confirms the outflow cavity nature of the latter, and from this point forward, we use the terms nebulosity and outflow cavity interchangeably.

The two cavities have almost exactly the same luminosity (see above) and similar sizes which suggests that the outflow is close to being in the plane of the sky (i.e. θ_incl ≳ 75–80°, with θ_incl defined from the line of sight) (e.g. Zhang & Tan 2018). The N1 cavity is mostly associated with red-shifted outflowing gas and is therefore pointing slightly away from us. 0n the other hand, N2 is mostly blue-shifted and is thus pointing slightly towards us. As the opening angle of the cavities is relatively large, see Sect. 5.3 for more details, it would be expected that both lobes contain a blue- and redshifted wing due to projection effects (Cabrit & Bertout 1986; Avison et al. 2021). Although weak, Fig. 7 shows hints of potential blue and red-shifted wings in both lobes even more than a beamsize away from the H1 fragment. Finally, the length of the ¹²CO(2–1) outflow, measured starting from the H1 fragment, is 0.66 pc for the red-shifted lobe and 0.55 pc for the blue-shifted lobe. These values are not corrected for inclination, but given the large inclination angle of the outflow these values are likely to be close to the actual lengths.

Although being much weaker than the one associated to the H1 fragment, a second identified bipolar outflow is located around the weak 70 μm source associated to H4 in the north-east of the map (see Fig. 7). No outflow has been detected towards the H2 and H3 fragments with the ¹²CO(2–1) data.

Table 2

Extracted fluxes in Jy for the H1 fragment and the two nebulosities from the Herschel, WISE and IRAC photometric maps.

Fig. 6

Spectral energy distributions obtained from Herschel, WISE and IRAC images, for the two nebulosities N1 and N2, and the H1 fragment.

Fig. 7

Top: integrated intensity contours of the redshifted (in red) and blueshifted (in blue) ¹²CO(2–1) high-velocity wings overlaid on the 70 μm map. The contours start at 10 K km s⁻¹ with increment of 2 K kms⁻¹. The green line and crosses indicate the locations of a spectra presented in the bottom panel with position 3 closely coinciding with the location of the H1 fragment. The white circle at the top indicates the beam size of the ¹²CO(2–1) data. Bottom: ¹²CO(2–1) spectra at the indicated locations in the top panel.

Table 3

Outflow momentum rate obtained for the blue- and redshifted parts of the two outflows.

5.2 The Outflow Momentum Rate

Integrating the ¹²CO(2–1) wing emission over the entire extent of the red- and blue-shifted outflows, one can estimate their C0 outflow momentum rate assuming the driving source is located at a velocity of −21 km s⁻¹. To calculate this, we used the technique described in Bontemps et al. (1996) and Duarte-Cabral et al. (2013) (see Appendix D for more details). The results from the calculations are listed in Table 3 for both the blue- and red-shifted components of the two bipolar protostellar outflows. To estimate the total outflow momentum rate, it is assumed that both sides eject an equal amount of momentum. The best estimate for the total outflow momentum rate is then obtained by taking two times the maximum value of the two sides (Duarte-Cabral et al. 2013). Assuming a mean inclination angle of 57.6° (Cabrit & Bertout 1992; Bontemps et al. 1996), the resulting momentum rate of the protostellar outflow originating in H1 is 2.8 × 10⁻³ M_⊙ km s⁻¹ yr⁻¹, which is similar to other CO estimates of the momentum rate at early stages of massive star formation (e.g. Duarte-Cabral et al. 2013; Tan et al. 2016).

Using $3 σ_{C^{18} O}$ $3{\sigma _{{{\rm{C}}^{18}}{\rm{O}}}}$ or $4 σ_{C^{18} O}$ $4{\sigma _{{{\rm{C}}^{18}}{\rm{O}}}}$ (instead of $3.5 σ_{C^{18} O}$ $3.5{\sigma _{{{\rm{C}}^{18}}{\rm{O}}}}$ ) when defining the outflow wing velocity ranges only changes the outflow momentum rate by +4% and -6%, respectively. This error is small compared to other sources of uncertainty in the determination of the momentum rate, such as: the inclination angle, the opacity of the C0 wings, the symmetry assumption, and the distance. As the axis of the protostellar outflow from H1 is expected to be located close to the plane of the sky, the assumed mean inclination could significantly underestimate the momentum rate. Working with a larger inclination of 70° or 80° increases the outflow momentum rate to 7.7 × 10⁻³ and 3.1 × 10⁻² M_⊙ km s⁻¹ yr⁻¹, respectively (see Table 3).

Finally, even though we have detected only one bipolar outflow centred on the H1 fragment, it is possible that multiple, unresolved, outflows are responsible for the observed outflow cavity. Only higher angular resolution observations of G345.88-1.10 will allow us to conclude on this particular question of outflow multiplicity.

5.3 Outflow Opening Angle

The opening angle of an outflow is often considered an indicator of how evolved the protostellar system is (e.g Arce et al. 2007; Kuiper et al. 2016), even though recent Hubble observations have questioned this idea (Habel et al. 2021). To determine the opening angle of the outflow emanating from the H1 fragment, we used both the ¹²CO(2–1) emission and the cavities seen with the 70 μm emission. The estimate of the opening angle was performed by creating a triangle connecting the outer parts of the outflow, which is demonstrated in Fig. E.1. This method assumes that the cavity is axi-symmetric. Measurements based on the ¹²CO(2–1) outflow give an estimated opening angle of 105° for the red-shifted outflow (N1) and 72° for the blue-shifted outflow (N2). Measurements based on the 70 μm emission give an estimated opening angle of 98° for N1 and 76° for N2. Despite the limited angular resolution of our ¹²CO(2–1) observations, the agreement between the two sets of measurements gives us confidence that the opening angle of the H1 outflow is about 90° ± 15°. These values are significantly larger than what is traditionally considered to be the upper limit of outflow opening angles for low-mass class 0 (i.e.~55°) and class I (i.e. ~75°) protostellar outflows (Arce & Sargent 2006). However, recent work by Habel et al. (2021) did present class 0 protostellar outflows with opening angles up to ~80°.

6 H_II Emission in the Cavities

In the previous section we showed that the presence of mid-infrared bright cavities at the centre of the G345.88-1.10 hub could be the consequence of protostellar outflow activity. In this section we investigate whether the G345.88-1.10 cavities could also be explained by the presence of ionising radiation from an embedded or undetected massive star (cluster). For that purpose, tracers of H II regions were used. Only two recent radio continuum surveys cover G345.88-1.10. The Parkes-MIT-NRAO (PMN) radio continuum survey (Griffith & Wright 1993) provides radio continuum data at 4.85 GHz (6.19 cm) at an angular resolution of 5′ (Condon et al. 1993). In this survey, G345.88-1.10 remains undetected with an upper limit of 100 mJy. This upper limit can be compared with the flux in the PMN survey for the 16 identified bipolar H_II regions near the Galactic plane between galactic longitudes ±60° (Samal et al. 2018; Deharveng et al. 2015). All these bipolar H_II regions but one have fluxes that are larger than the flux upper limit for G345.88-1.10 by at least one order of magnitude. 0nly G051.61-0.36 from Samal et al. (2018) has a flux of ~0.2 Jy, which is still a factor two above the upper limit on G345.88-1.10. Having a closer look at the 70 μm emission of G051.61-0.36 in the HiGAL maps (Molinari et al. 2010), this source looks nothing like G345.88-1.10 as it is a 70 μm point source without 70 μm emission in the bipolar cavities. A comparison with this set of bipolar Hii regions thus suggests that G345.88-1.10 is not a classical bipolar H II region. The other radio continuum survey that covers G345.88-1.10 is the Sydney University Molonglo Sky Survey (SUMSS) at 843 MHz with an angular resolution of 45″ (Bock et al. 1999). In this survey some very weak radio continuum emission is detected towards G345.88-1.10. (see Fig. 8). This emission is not centred on the massive H1 fragment but is observed towards both cavities. This emission is very weak, around three times the noise level, which explains why G345.88-1.10 is missing in the SUMSS source catalogue (Mauch et al. 2003). Measuring the SUMSS flux of both cavities we obtain 26 mJy for N1 and 21 mJy for N2. However, it should be noted that G345.88-1.10 is located ~30′ away from a very bright SUMSS source. As a result, a complex low-intensity background of ~10 mJy beam⁻¹ can be seen in the vicinity of G345.88-1.10, contributing to an overestimate of the cavity fluxes presented above.

Lastly, we inspect the H30α recombination line covered by the APEX setup. In order to improve the signal-to-noise ratio we resampled the data to 2 km s⁻¹ velocity channels. However, no emission is detected anywhere in the region (see Appendix F). Averaging the spectra over both cavities gives a 3σ upper-limit for the line of 63 and 110 mK for the N1 and N2 lobes, respectively. Assuming a typical electron temperature of 7500 K at the galactocentric radius of G345.88-1.10 (Quireza et al. 2006), the equations presented in Appendix F allow us to calculate an upper limit for the cavities’ electron density. For N1 and N2 we obtain n_e < 3.4 × 10² cm⁻³ and n_e < 4.6 × 10² cm⁻³, respectively. Combined with the assumed electron temperature, these electron density upper limits translate into thermal pressure upper limits for the ionised gas of 5.0 × 10⁶ and 6.9 × 10⁶ K cm⁻³ in N1 and N2, respectively. Compared to the thermal and ram pressure of the molecular gas, estimated to be ~2.6 × 10⁵ K cm⁻³ (assuming T_K = 20 K) and 1.2 × 10⁷ K cm⁻³ (assuming v = 2 km s¹) respectively, the thermal pressure upper limit of the ionised gas is below the ram pressure in the hub. This indicates that the clearing of gas and dust at the location of the nebulosities cannot be due to the ionising radiation of an unresolved stellar cluster. This hints at the idea that the protostellar outflow might play a role in the formation of the observed bipolar cavity centred on H1. Using the upper limit at 5 GHz, which typically is optically thin, with the equation from Schmiedeke et al. (2016) gives a more stringent n_e < 2.5 × 10² cm⁻³ and resulting thermal H_II pressure <3.7 × 10⁶ K cm⁻³.

The H30α line and 5 GHz radio continuum thus only provide an upper limit on the pressure from ionised gas, and the 843 MHz continuum emission cannot be used to estimate the electron density, even though it is detected, because the emission can be optically thick at this wavelength. Furthermore, at these long wavelengths the emission might even be the result of synchrotron radiation from a powerful jet associated with the outflow (e.g. Carrasco-González et al. 2010; Rodriguez-Kamenetzky et al. 2017). Dedicated radio continuum observations are required to resolve this question.

Fig. 8

Radio continuum image of G345.88-1.10 at 843 MHz from the SUMSS survey (Bock et al. 1999). The 70 μm emission observed towards the bipolar cavity is highlighted with black contours (2000, 4000, 6000, 8000 and 10 000 MJy sr⁻¹).

7 Radiative Transfer Calculations

The mid-infrared images of G345.88-1.10, presented in Figs. 3 and 5, show at its centre a large bright bipolar nebulosity which is closely associated to a bipolar outflow emanating from the H1 fragment. Infrared bright outflow cavities are commonly observed in the near-infrared, typically at a wavelength of a couple of microns, where photons are efficiently scattered by dust grains located on the surface of the cavity. In the case of G345.88-1.10, the cavity is the brightest at 70 μm, a wavelength at which scattering is negligible. Therefore, the mid-infrared dust emission towards the outflow cavities of G345.88-1.10 must be thermal dust emission. The obvious mechanism to heat up dust grains on the walls of an outflow cavity is by direct illumination of a central protostar. The most probable location of that protostar is embedded within the infrared dark H1 fragment.

To investigate whether radiative heating from a deeply embedded source could explain the observed mid- to far-infrared emission of G345.88-1.10, we performed a series of dust continuum radiative transfer calculations.

7.1 Pandora Setup

In this paper we use the python-based radiative transfer framework Pandora (Schmiedeke et al. 2016). Here we briefly summarise the components embedded in Pandora, that we employ in our study. We use RADMC-3D (Dullemond 2012) for a self-consistent determination of the dust temperature and calculation of synthetic continuum maps, and MIRIAD (Sault et al. 1995) for the post-processing. We explore the parameter space by hand.

7.1.1 Mesh refinement of RADMC-3D

For the radiative transfer calculations, we use a square box of 6 pc in size containing a dust distribution and heating sources. In order to resolve dust over-densities present within the box volume, adaptive mesh refinement (AMR) is used (Berger & Oliger 1984; Berger & Colella 1989) by RADMC-3D. Here, the AMR grid starts with 11 cells along every major axis and is refined into 8 subcells at each refinement level. A cell is refined if the density difference within a cell exceeds 10%. Then, a minimum cell size is specified, which, in our case, is set to 100 AU to properly resolve all dust distributions that will be studied. For each grid the dust temperature, and therefore the thermal dust emission, is then self-consistently computed using the Monte Carlo method of Bjorkman & Wood (2001).

7.1.2 Dust Continuum Images

When the photon packages leave the volume grid, they produce an image at different wavelengths in units of Jy pixel⁻¹. For each of the simulations 10⁷ photon packages were used (Schmiedeke et al. 2016). To compare these synthetic images with the observations of G345.88-1.10, Miriad is employed (Sault et al. 1995) in Pandora. Using the beam size at every wavelength, it convolves the artificial image with a Gaussian beam. This returns the convolved image in Jy beam⁻¹ which can be easily converted to MJy sr⁻¹ units and compared with the observations.

7.2 Model Setup

7.2.1 Heating Sources

The central heating source (i.e. the embedded protostar) is modelled as a blackbody with a specific luminosity. Three values were used for this: 500, 2000 or 4000 L_⊙. In addition to this internal source of heating, the RADMC-3D calculations include an interstellar radiation field (ISRF) which consists of the cosmic microwave background (CMB), the far infrared emission from dust grains (Draine & Li 2007), the starlight background from Draine (2011), and a UV background from Mathis et al. (1983).

7.2.2 Dust Density Distribution

Because of the uncertainties related to the unresolved morphology of the H1 fragment, we here used combinations of three different types of density distributions.

First, we used a modified Plummer model as described in Qin et al. (2011) $n (r) = \frac{n_{c}}{{(1 + {| r |}^{2})}^{η / 2}}$ $n\left( {\bf{r}} \right) = {{{n_c}} \over {{{\left( {1 + {{\left| {\bf{r}} \right|}^2}} \right)}^{{\eta \mathord{\left/ {\vphantom {\eta 2}} \right. \kern-\nulldelimiterspace} 2}}}}}$ (3)

where n_c is the central density. The power law index η was taker to be 2 (Schmiedeke et al. 2016), and |r| is given by $| r | = \sqrt{{(\frac{x}{r_{0, x}})}^{2} + {(\frac{y}{r_{0, y}})}^{2} + {(\frac{z}{r_{0, z}})}^{2}}$ $\left| {\bf{r}} \right| = \sqrt {{{\left( {{x \over {{r_{0,x}}}}} \right)}^2} + {{\left( {{y \over {{r_{0,y}}}}} \right)}^2} + {{\left( {{z \over {{r_{0,z}}}}} \right)}^2}}$ (4)

where r_0,x, r_0.y, r_0,z are the characteristic sizes of the Plummer model in the direction of the x, y and z axes. The Plummer model is specifically used to model the core or cloud embedding the protostar.

The second density distribution is a passive flared disk (PFD) which follows the description in Pineda et al. (2011) $n (R, z) = n_{0} {(1 + \frac{r_{0}}{R})}^{β - q} \exp (- 0.5 {(\frac{z}{h (R)})}^{2})$ $n\left( {R,z} \right) = {n_0}{\left( {1 + {{{r_0}} \over R}} \right)^{\beta - q}}\exp \left( { - 0.5{{\left( {{z \over {h\left( R \right)}}} \right)}^2}} \right)$ (5)

with: $h (R) = h_{0} {(\frac{R}{r_{0}})}^{β}$ $h\left( R \right) = {h_0}{\left( {{R \over {{r_0}}}} \right)^\beta }$ (6)

In these equations, R is the radial distance in the equatorial plane, n₀ is the characteristic density, r₀ the characteristic radius, β the disk flaring power, q the surface density radial exponent, and h₀ the scale height at the characteristic radius. This PFD starts at an inner radius r_i.

Lastly we also take the presence of outflow cavities into account, which are characterised by a length, an opening angle and a constant density. All models use a cavity opening angle of 80° and a length of 10⁵ AU (-0.5 pc).

Combinations of these three density distributions were used to produce the complete density distribution of the radiative transfer models. The models will be specified in more detail in the next section. When combining the different density distributions in a single model for G345.88-1.10, the total density in each cell ‘j’ was determined by $n_{j} = \sum_{i = 1}^{N} n_{i, j} (r)$ ${n_j} = \sum\limits_{i = 1}^N {{n_{i,j}}({\bf{r}})}$ (7)

where N is the number of dust distributions used for the model, and n_i,j the density of each density distribution ‘i’ in the cell ‘j’. This expression is not valid for the outflow cavity regions, where only the constant density of the cavity contributes.

Table 4

Density distribution elements for all three models.

7.3 Model Grid and Parameter Space

To study the observed features in the Herschel continuum images of G345.88-1.10 around the H1 fragment, a set of potentially representative radiative transfer models were run. As mentioned before, this set of models uses a central heating source that is embedded in a dust distribution at the origin of the x, y and z axes. The dust distributions for this set of models can be divided into 3 categories, and are from here on referred to as Model 1, 2 and 3. Their composition is presented in Table 4 and Fig. 9.

For each of the three models, we probed a range of values for most of the parameters in Eqs. (3)–(6). These are summarised in Table 5. We performed three simulations for each parameter combination, one for a central luminosity of 500 L_⊙, 2000 L_⊙ and 4000 L_⊙ each, leading to a total of 513 different RADMC-3D radiative transfer calculations: 243 setups were run for Model 1, 108 setups for Model 2, and 162 setups for Model 3. The radiative transfer calculations thus cover a wide range of possible configurations. As a consequence of the expected high inclination angle for the G345.88-1.10 outflow cavity with respect to the line of sight, an inclination of 90° was given to all the models. When running models where the cavities have an inclination of 80°, it is found, as expected (Johnston et al. 2011; Zhang & Tan 2018), that there is a substantial difference in brightness for the two cavities. This again indicates that the cavities in G345.88-1.10 have a very small angle with respect to the plane of the sky. Note that we limited ourselves to a maximum luminosity of 4000 L_⊙ for several reasons: The measured bolometric luminosity of the entire system (cavities and H1 fragment) is about 4000 L_⊙ (see Sect. 4.2); above a luminosity of 4000 L_⊙ increasingly bright central H_II regions are expected to form (Kalcheva et al. 2018; Csengeri et al. 2018) while we do not detect any in G345.88-1.10; models with L = 4000 L_⊙ are found to be already too luminous in order to reproduce some key features (see next section).

All calculations make use of the thin ice mantle dust grain model as there evidently is some heating in this region (Ossenkopf & Henning 1994). But it is likely that, given the range of densities and temperatures within the system, one dust grain model is not enough. This is one of the main limitations of the work we present here.

Fig. 9

Schematic representation of the three density models used for the radiative transfer calculations. The blue ellipse/circle represents the Plummer (flattened) dense core in Model 1 & 3, and the Plummer cloud for Model 2. The grey triangle represent the cavities, while the green ring represents the dense central disk. Finally, the central yellow star symbol represents the embedded heating source.

Table 5

Input parameters that were used for the grid of radiative transfer simulations (Bontemps et al. 2010; Palau et al. 2013; Beuther et al. 2015a; Beltrán & de Wit 2016), subdivided based on the three different model setups.

7.4 Results

7.4.1 The 70 μm Bright Cavities

Since the unusually bright 70 μm cavities is one of the most peculiar aspects of G345.88-1.10, we first focus on whether or not it is possible to produce a bipolar nebulosity that is brighter at 70 μm than the central core itself. Across all 513 radiative transfer calculations we made, about 12% have their 70 μm brightness peaking within the cavities, highlighting the fact that only a small subset of radiative transfer calculations manage to reproduce this unusual feature. The number of runs that do manage to produce a dominant 70 μm nebulosity is non-uniform across the different models. Model 2 and 3 are more successful at it (see Table 6) which shows that a dense compact object, a flared disc in our case, is required to produce 70 μm nebulosities through radiative heating. It can also be inferred from Table 6 that, within the studied range, the occurrence of nebulosities has little dependence on the source luminosity.

Analysing the runs that do produce a 70 μm-bright nebulosity in more detail, we realise that none of them manages to match the observed brightness of the G345.88-1.10 nebulosity. The run with the brightest nebulosity only recovers about 15% of the observed 70 μm intensity for G345.88-1.10. Furthermore, more than 80% of the nebulosities in the RT calculations produce less than 2.5% of the observed flux in G345.88-1.10. These low fractions suggest that the observed nebulosities around H1 cannot be solely the result of radiative heating.

The 70 μm dust continuum image with the brightest modelled nebulosity is shown in Fig. 10, and was achieved by a run with Model 3. This image can be directly compared to Fig. 3 (left), the observed 70 μm image of the G345.88-1.10 nebulosity. From this comparison we first notice that there is an order of magnitude difference in brightness, G345.88-1.10 being the brighter of the two. We also find that, while the overall sizes are in relatively good agreement, the nebulosity morphologies are different. Most noticeably, the modelled nebulosities have the highest flux density at the opening of the outflow cavity where dust heating is the most efficient, while the observed nebulosities in G345.88-1.10 have a more complex distribution. This is another indication that radiative heating might not be the driving mechanism behind the observed cavities in G345.88-1.10, although the more complex cavity morphology of the region, compared to the models, likely also plays a role in this.

Finally, in G345.88-1.10 the cavity is still bright at 160 μm, see Fig. 5. The resulting 70-160 μm flux ratios of the G345.88-1.10 cavities are found to be 2.8 and 2.5 for N1 and N2, respectively (see Table 2). However, in all radiative transfer calculations, the modelled nebulosity is always brighter at 160 μm, with maximum 70–160 μm flux ratios of 0.45, 0.58, and 0.65 for the 500 L_⊙, 2000 L_⊙ and 4000 L_⊙ runs, respectively.

All these results show that there are important differences between the observed nebulosities in G345.88-1.10 and those obtained within our radiative transfer calculations. This implies that even though the flashlight effect, giving rise to an anisotropic thermal radiation field (Yorke & Sonnhalter 2002; Kuiper et al. 2010), could give rise to weak 70 μm nebulosities around young protostars, it cannot explain the far/mid-infrared properties of the outflow cavities that are observed towards G345.88-1.10 in the performed simulations.

Table 6

Fraction of radiative transfer 70 μm nebulosities appearing in a perfectly edge-on configuration (in %) as a function of luminosity and density distribution.

Fig. 10

70 μm map of the radiative transfer simulation that produced the brightest nebulosities. These cavities are produced by ‘Model 3’ with a central object luminosity of 4000 L_⊙, a characteristic radius of 10⁴ AU, a characteristic density of 10⁷ cm⁻³ and a cavity density of 10⁴ cm⁻³. This is to be directly compared to Fig. 3 (left). Note the order of magnitude difference in flux density between the two images.

7.4.2 The 70 μm quiet fragment

In Fig. 5, the H1 fragment is best identified from 250 μm upwards. At shorter wavelengths, dust continuum emission traces either both the cavities and the core (160 μm) or the cavities only (from 70 μm downwards). In order to test the ability of the radiative transfer calculations to reproduce the emission properties of the H1 fragment, we thus focus on the wavelength range between 250 and 500 μm which also includes the Herschel 350 μm data. The goodness of a specific run is estimated via a reduced χ² calculation: $χ_{core}^{2} = \frac{1}{N} \sum_{λ} \frac{{(S_{λ}^{core} - O_{λ}^{core})}^{2}}{σ_{λ}^{2}}$ $\chi _{{\rm{core}}}^2 = {1 \over N}\sum\limits_\lambda {{{{{\left( {S_\lambda ^{{\rm{core}}} - O_\lambda ^{{\rm{core}}}} \right)}^2}} \over {\sigma _\lambda ^2}}}$ (8)

where $S_{λ}^{core}$ $S_\lambda ^{{\rm{core}}}$ are the core fluxes estimated from the synthetic images at the wavelength λ (i.e. 250, 350 and 500 μm, respectively). $O_{λ}^{core}$ $O_\lambda ^{{\rm{core}}}$ are the observed fluxes towards the H1 fragment (see Table 2). N is the number of wavelengths (i.e. 3), and σ_λ is the uncertainty on the observed flux at the wavelength λ (i.e. 20% of the observed flux), to take into account the distance and calibration uncertainties.

Figure 11 shows the corresponding $χ_{core}^{2}$ $\chi _{{\mathop{\rm core}\nolimits} }^2$ values as a function of the central source luminosity. Only models achieving $χ_{core}^{2} < 3$ $\chi _{{\mathop{\rm core}\nolimits} }^2 < 3$ are displayed. At first sight, there appears to be no clearly preferred luminosity, both 500 L_⊙ and 4000 L_⊙ runs achieve low $χ_{core}^{2}$ $\chi _{{\mathop{\rm core}\nolimits} }^2$ values. However, at 4000 L_⊙ these low $χ_{core}^{2}$ $\chi _{{\mathop{\rm core}\nolimits} }^2$ values are only obtained for three specific setups of model 2 (which contains only a disk and ambient cloud around the central heating source), while at 500 L_⊙ low values are obtained for all three models. Interestingly, we can see that while Model 1 and 3 runs are doing a better job at low luminosity and a worse job at high luminosity, Model 2 runs do not show this tendency. This is a direct consequence of the reprocessing of the radiation from the central protostar by dust. In Model 2, photons can escape the system almost freely and therefore contribute significantly less to heat up the dust within the disc, leading to far-infrared fluxes that are compatible with an infrared quiet core. However, for Models 3 and 1, the dense central Plummer core absorbs and re-emits much of the source radiation. At high luminosity this translates into a central core that is too bright compared to the observed fluxes. As a result, Models 3 and 1 runs only manage to reproduce the 250, 350 and 500 μm emission properties of H1 at low luminosity.

At 70 μm, the complexity of the nebulosity emission makes it hard to determine a 70 μm flux for the H1 fragment. However, we can still derive an upper limit for the 70 μm flux of the H1 fragment which is 5 Jy. We can now check whether the radiative transfer calculations that provide the best fits to the emission properties of the core at long wavelengths are also compatible with this 70 μm flux upper limit. To select the radiative transfer calculations that fit reasonably well the observations of the H1 fragment in the 250–500 μm range, we consider models for which $χ_{core}^{2} \leq 2$ $\chi _{{\mathop{\rm core}\nolimits} }^2 \le 2$ . This provides a set of 54 simulations to work with, 37 of which have a central luminosity of 500 L_⊙, 9 have a central luminosity of 2000 L_⊙ and 8 have a central luminosity of 4000 L_⊙. Within this subset of models, the 2000 L_⊙ and 4000 L_⊙ runs have a minimum 70 μm flux of 81 Jy and 29 Jy, respectively. This is more than a factor of 6 higher than the observed upper limit. For the 500 L_⊙ runs, the distribution of obtained 70 μm fluxes are plotted in Fig. 12. These radiative transfer calculations produce 70 μm fluxes that are low enough to be compatible with the observed flux upper limit. This shows that the central protostellar source embedded within H1 must have a luminosity ≤500 L_⊙. However, for such low luminosity models, the discrepancy at 70 μm between the observed and modelled cavities increases to several orders of magnitude (see Appendix G). This implies that our radiative transfer calculations cannot reconcile the presence of an infrared quiet core with a luminous mid-infrared bright cavity.

Lastly, it should be noted that we used the same thin ice mantle dust grain model for all calculations. Changing the dust distribution to include larger grain mantles can increase the emissivity up to a factor three at 70 μm (Ossenkopf & Henning 1994). However, this is insufficient to resolve the orders of magnitude discrepancy with the observations of G345.88-1.10.

Fig. 11

$χ_{core}^{2}$ $\chi _{{\mathop{\rm core}\nolimits} }^2$ as a function of the central source luminosity for all models with $χ_{core}^{2} < 3$ $\chi _{{\mathop{\rm core}\nolimits} }^2 < 3$ . Model 1: red circle; Model 2: blue triangle; Model 3: green square.

Fig. 12

70 μm flux distributions for the simulated cores with $χ_{core}^{2} \leq 2$ $\chi _{{\mathop{\rm core}\nolimits} }^2 \le 2$ for an internal luminosity of 500 L_⊙. The vertical dashed line indicates the upper limit of 5 Jy observed towards G345.88-1.10.

8 Discussion

In the following discussion, we consider that the major axis of the G345.88-1.10 bipolar cavity and its associated outflow have a large inclination angle with respect to the line of sight. As discussed in previous sections, this is mostly supported by the similar morphology and brightness of the two lobes of the nebulosity.

8.1 A Potentially very Young Protostar Generating Powerful Mass Ejection

All four fragments identified within G345.88-1.10 are massive and quiet in the mid-infrared. This includes the H1 fragment for which radiative transfer calculations show that the luminosity of the embedded protostar has to be ≤ 500 L_⊙. The same radiative transfer calculations also show that radiative heating from such a low luminosity source appears to be unable to account for the order of magnitude brighter cavities. Additionally, with the H30α recombination line and the 5 GHz radio continuum we do not detect a thermal H II region. Therefore, we consider here the possibility that the heating in the cavity is mostly mechanical, originating from shocks of a powerful protostellar outflow. In the following we investigate the plausibility of such a scenario.

In the case where the cavity brightness is powered by the outflow kinetic energy, then the cavity luminosity could provide an estimate of the outflow mechanical power. This is however complicated by the fact that shock heating is generally not radiated away by dust continuum emission. It could be considered that emission from shock cooling lines like [0I] contributes to the total observed continuum emission (e.g. in the Herschel 70 μm band). This is however very unlikely as this would require excessively bright cooling lines. Recently developed shock models did however show that up to 30% of shock kinetic flux can escape as FUV Lyα photons (Lehmann et al. 2020, 2022). These photons could then heat the dust and excite cooling lines within the cavities.

In a scenario where the observed cavity luminosity is powered by the mechanical heating of the protostellar outflow, one may link luminosity and outflow momentum rate via $L = \frac{1}{2} \frac{f_{irr}}{f_{ent}} F_{{CO}^{v} jet}$ $L = {1 \over 2}{{{f_{{\rm{irr}}}}} \over {{f_{{\rm{ent}}}}}}{F_{{\rm{C}}{{\rm{O}}^v}{\rm{jet}}}}$ (9)

where v_jet is the jet velocity. To correct for the fact that only a fraction of the outflow mechanical energy produces Lyα photons (Lehmann et al. 2020, 2022) that can heat the dust, we include the factor f_irr. Here, we assume that 20% of the outflow energy is being radiated away. Further, we use a C0 entrainment efficiency of 50% (i.e. f_ent; Duarte-Cabral et al. 2013) and the observed C0 momentum rate at an 80° inclination angle (F_CO). This gives a jet velocity of 3.8 × 10³ km s⁻¹. This is significantly above the typical velocity range of 100–1000 km s⁻¹ that is reported for protostellar jets (e.g. Torrelles et al. 2011; Gusdorf et al. 2017; Anglada et al. 2018). It thus appears that mechanical heating by a protostellar outflow does not provide a straightforward explanation for the high luminosity of the cavities.

One should keep in mind that the different parameters (i.e. distance, inclination angle, and CO entrainment efficiency) that enter the jet velocity estimate are highly uncertain. Towards some high-mass star forming regions, it was found that part of the outflow is also detected in [C_II] (e.g. Schneider et al. 2018; Sandell et al. 2020; Bonne et al. 2022). This implies that a fraction of the outflow mechanical power might be completely missed by CO observations. If the CO entrainment efficiency is lower than 0.5, this could reduce the required jet velocity towards a more compatible value compared to other protostellar jet velocities. This argument is particularly important if Lyα emission is produced by the shocks as these photons can dissociate CO.

Finally, if mechanical heating were to be the origin of the source’s luminosity, then one might expect evidence of ongoing shocks on the cavity walls. Unfortunately, we do not have shock tracer observations at the moment, but the C¹⁸O(2–1) and ¹³CO(2–1) spectra observed towards G345.88-1.10 present several skewed line shapes (see Fig. 13 and Appendix C). A closer inspection of these skewed spectra shows that they are not randomly positioned, but instead are prominent in, and at the edges of, the two lobes of the bipolar nebulosity, suggesting a direct link between skewed line profiles and the cavities. These skewed spectra, especially ¹³CO(2–1), show a notable resemblance in shape to the spectra observed with the shock tracer SiO(2–1) around IR quiet massive dense cores in Cygnus-X with both a narrow and broader components (Motte et al. 2007; Duarte-Cabral et al. 2014). In G345.88-1.10, the velocity range for C¹⁸O(2–1) and ¹³CO(2–1) cover up to 8–10 km s⁻¹. This is generally smaller than in Motte et al. (2007) where the tail covers a velocity interval up to 20-30 km s^-1. However, such a difference could easily be explained by the high inclination angle of the outflow axis with respect to the line-of-sight. Even though the shape of the ¹³CO(2–1) and C¹⁸O(2–1) spectra towards the cavity are suggestive of the presence of shocked gas, it should be noted that with current observations it is not possible to exclude that these spectra might also be associated with strong line-of-sight velocity gradients. High-angular resolution observations of well-known shock tracers will be necessary to resolve this issue.

Fig. 13

¹²CO(2–1), ¹³CO(2–1) and C¹⁸O(2–1) spectra extracted in the N1 (top) and N2 (bottom) cavities. It is observed that the high-velocity outflow wing, seen in ¹²CO(2–1), are at the opposite side of the spectrum with respect to the lower-velocity tails that are also seen in ¹³CO(2–1) and C¹⁸O(2–1) on top of ¹²CO(2–1).

8.2 Global Collapse of a Hub

Mass rates from infall and protostellar outflows are connected to each other. While in the low-mass case, it is believed that proto-stellar outflows are powered by the collapse of protostellar cores, in the high-mass case there are indications that it is the collapse of the entire parsec-scale parent clump (e.g. Avison et al. 2021). In a scenario where the observed mid-infrared cavities might be shaped by the H1 outflow, a link between clump collapse and outflows suggests a link between clump collapse and the presence of cavities. In this section, we investigate the kinematics of the G345.88-1.10 hub in an attempt to characterise its collapse.

Figure 2 shows the ¹²CO(2–1), ¹³CO(2–1), and C¹⁸O(2–1) spectra averaged over the extent of the G345.88-1.10 hub. While the former two spectra show a double-peaked spectrum with predominant blue-shifted spectral emission, the latter shows a single-peaked line centred near the dip of the other two. This is the well-known blueshifted self-absorbed spectral line shape which is often interpreted as a signature of collapse (e.g. Mardones et al. 1997). The CO isotopologues are not the most commonly used tracer of gravitational inflow as they tend to trace more external layers of the clouds. However, in cases where clouds are self-gravitating and collapsing on large scales, these tracers can be used to probe infall velocities (Schneider et al. 2015b). Using the simple analytical model of Myers et al. (1996), a first estimate of the line-of-sight infall velocity can be made from the ¹³CO(2–1) spectrum, which gives a derived infall velocity of 0.5 km s⁻¹ (using ¹²CO(2–1) gives a similar value of 0.6 km s⁻¹). However, De Vries & Myers (2005) noted that Myers’ model tends to underestimate the infall velocity by a factor 2, leading to a corrected infall velocity v_inf ~ 1 km s⁻¹. As a result of the large opacity of CO lines and the depletion of CO isotopologues at moderate densities (e.g. Tafalla et al. 2004), the derived infall velocity is most likely representative of the dynamics of the hub outer layers only.

Figure 14 presents a map of the C¹⁸O(2–1) hub centroid velocity obtained in C¹⁸O(2–1). In that plot it can be clearly seen that velocity gradients show a rapid increase around the H1 and H2 fragments, in the densest region of the hub. This behaviour has been proposed to be the result of gravitational collapse by several groups (e.g. Hartmann & Burkert 2007; Peretto et al. 2007, 2013; Smith et al. 2013; Gόmez & Vázquez-Semadeni 2014; Hacar et al. 2017; Watkins et al. 2019). The velocity gradients near the H1 and H2 fragments range from 4 to 6 km s⁻¹ pc⁻¹, while the average value over the hub is ${1.9}_{- 0.2}^{+ 0.3}$ $1.9_{ - 0.2}^{ + 0.3}$ km s⁻¹ pc⁻¹. Velocity gradients of this magnitude and larger have been observed at significantly smaller scales (e.g. Beuther et al. 2015b; Dhabal et al. 2018). However, at the spatial resolution of the APEX data (i.e. ~ 0.3 pc) these gradients are large when comparing to other star forming regions, whether they are low-mass or high-mass (e.g. Schneider et al. 2010; Henshaw et al. 2013, Henshaw et al. 2014; Peretto et al. 2014; Bonne et al. 2020). This indicates that the G345.88-1.10 hub is particularly dynamic overall, and gas might be accelerated towards its centre, as expected in the case of the global collapse of centrally condensed clouds (e.g. Peretto et al. 2007). However, another origin for organised gas motion is rotation. We therefore checked whether the observed velocity field could be explained by rotation alone. But assuming both spherical and flattened geometries, the rotational energy for the hub is less than 15% of the hub gravitational energy. Thus, even though we cannot rule out its presence, it is clear that rotation on its own cannot prevent the collapse of the G345.88-1.10 hub. It could be that the hub velocity field is the result of both rotation and kinematics driven by a strong magnetic fields that provides support against collapse (e.g. Inoue et al. 2018; Bonne et al. 2020; Arzoumanian et al. 2022). However, a simpler explanation would be that gravity is the main driver of the observed velocity field, as suggested by the infall spectral signatures discussed above.

To further investigate the role of gravity one can compare the observationally-derived inflow timescale to the theoretical free-fall time. The inflow timescale is given by: $t_{\inf} = \frac{R}{v_{\inf}}$ ${t_{\inf }} = {R \over {{v_{\inf }}}}$ (10)

where R is the hub radius. Taking R = 1 pc and v_inf = 1 km s⁻¹, we obtain t_inf = 1 Myr. On the other hand, the free fall time scale is given by: $t_{ff} = \sqrt{\frac{3 π}{32 G ρ}}$ ${t_{{\rm{ff}}}} = \sqrt {{{3\pi } \over {32G\rho }}}$ (11)

where ρ is the average mass density within the hub. With a mass of 2900 M_⊙ within a radius of 1 pc (see Sect. 3.2), we obtain ρ = 4.6 × 10⁻²⁰ g cm⁻³ which corresponds to a molecule number density of $n_{H_{2}} = 1.2 \times 10^{4}$ ${n_{{{\rm{H}}_2}}} = 1.2 \times {10^4}$ cm⁻³. Plugging this density in the equation of the free-fall time we get t_ff = 0.3 Myr. We can conclude that the collapse time of the G345.88-1.10 hub is about 3 times longer than its free-fall time. This implies that the collapse of G345.88-1.10 is not pressure-free, but it is still significantly faster than what we expect for the quasi-static evolution of clumps (e.g. Krumholz & Tan 2007). Now, assuming that the mass inflow is spherically symmetric we can compute the hub mass infall rate using: ${\dot{M}}_{\inf}^{hub} = 4 π R^{2} v_{\inf} ρ$ $\dot M_{\inf }^{{\rm{hub}}} = 4\pi {R^2}{v_{\inf }}\rho$ (12)

Here we use the same values as those used for the timescales derivation, except for the radius where we use R = 0.5 pc. The reason for this is that the average hub density p is reached at a radius that is about half the external hub radius of 1 pc. With these values in hand we obtain ${\dot{M}}_{\inf}^{hub} = 2.2 \times 10^{- 3} M_{⊙} {yr}^{- 1}$ $\dot M_{\inf }^{{\rm{hub}}} = 2.2 \times {10^{ - 3}}\,{M_ \odot }\,{\rm{y}}{{\rm{r}}^{ - 1}}$ . This is very similar to what has been derived for the prototypical SDC335 hub filament system (Peretto et al. 2013).

The gravitational collapse of the G345.88-1.10 hub provides the means for rapid mass provision to the centre, along with non-steady accretion flows which might open the cavity on short timescales through an accretion burst. This can be the result of e.g. jet precession (Rosen & Krumholz 2020) and opening of the jet itself (Cesaroni et al. 2018). The modelling by Cesaroni et al. (2018) indicates that the jet opening angle can increase to values as high as 40–50°, but it should also be noted that this is model dependent. The simulations by Rosen & Krumholz (2020) indicate that a jet precession can be very rapid at the early stages of high-mass star formation with rapid angle changes that on average reach 40° kyr⁻¹ with apparent peaks as high as 250° kyr⁻¹. However, it should be noted that the opening of the cavities by an accretion burst also faces the issue that this is generally associated with a strong increase of the core brightness in the far-infrared (e.g. Stecklum et al. 2021; Hunter et al. 2021), which is not seen in G345.88-1.10. It could be that, after an accretion burst has occurred, there is a time delay between the dimming of the core and the dimming of the cavity. Such a time delay would be of the order of the jet’s crossing time of the cavity, i.e. ≲1000 yr. We could therefore be catching G345.88-1.10 right at that moment where the cavity is still bright but the core has dimmed. Although this would be a noteworthy coincidence, it can fit with the cooling timescale we estimate for the massive H1 fragment which should be lower than 20 yr and can be as low as a couple of days or months, see Appendix H. These timescales also seem to be in agreement with surveys of continuum variability, at least in low-mass star forming cores that typically find very rapid changes on the timescales of weeks to several years (e.g. Johnstone et al. 2013; Mairs et al. 2017; Lee et al. 2021). Because of these very short cooling timescales, it is even possible that multiple accretion bursts spread over few 10-100 yr are at the origin of the heated cavities. In Appendix H, we also estimated the cooling timescale in the cavities for which we find values between 1.7 × 10³ and 2.7 × 10⁴ yr which is significantly longer than for the core. Furthermore, the lower end of this estimated timescale is of the same order of magnitude as the crossing timescale for the jet in the cavities. This might explain why the full cavity is bright instead of only several subregions in the cavity.

Fig. 14

Velocity centroid map obtained by fitting a Gaussian profile to the C¹⁸O(2–1) spectra. The calculated velocity gradient is overplotted on top of it, using streamplot (https://matplotlib.org/stable/api_as_gen/matplotlib.pyplot.streamplot.html), showing organised motion. The colourbar of this streamplot, located on top of the figure, indicates the magnitude of the velocity gradient at each location. The black circle indicates the beam size of the presented APEX observations and the black crosses indicate the position of the 4 massive fragments.

9 Summary and Future Prospects

We have presented the first study of the G345.88-1.10 hub filament system ( $d = {2.26}_{- 0.21}^{+ 0.30}$ $d = 2.26_{ - 0.21}^{ + 0.30}$ kpc, M = 2.9 × 10³ M_⊙ within a radius of 1 pc). The analysis in this paper focuses on the massive 70 μm-quiet H1 fragment, located at the hub centre, and the associated bright bipolar cavity. We combined archival Herschel, WISE, Spitzer and 2MASS data with APEX observations of ¹²CO(2–1), ¹³CO(2–1), C¹⁸O(2–1) and H30α towards the hub. The main results are:

The H1 fragment has an estimated mass of ${2.1}_{- 1.1}^{+ .2} 10^{2}$ $2.1_{ - 1.1}^{ + - .2}{10^2}$ M_⊙ within a radius R_eff = 0.14 pc, and a luminosity upper limit of 10² L_⊙ from direct integration of the SED.
Bipolar infrared nebulosities are centred on the H1 fragment. The two lobes of the nebulosity have similar overall shapes and luminosities, with a combined luminosity ~4.4 × 10³ L_⊙. It is believed that the main axis of the bipolar nebulosity is lying close to the plane of the sky.
¹²CO(2–1) high velocity wings, centred on the H1 fragment, are observed in the bipolar nebulosities with an estimated momentum rate of 3.1 × 10⁻² M_⊙ km s⁻¹ yr⁻¹. The nebulosities thus trace the outflow cavity of a powerful protostellar outflow emanating from the H1 fragment with an opening angle of ~90 ± 15°.
Weak 843 MHz radio continuum emission is detected towards the cavities in the SUMSS survey. However, H30α and 5 GHz radio continuum emission is not detected. This puts upper limits on the thermal pressure of potential H_ii gas. As a result, ionising radiation appears unable to open the cavities.
Radiative transfer calculations of a disc/core/cavity system confirms that the luminosity of the embedded source in the massive H1 fragment has to be ≤500 L_⊙. Most importantly, these calculations also show that the bright 70 μm nebulosities at the centre of G345.88-1.10 cannot be explained by the radiative heating from an embedded object.
Potential heating of the cavities by shocks from the powerful protostellar outflow also appears unable to explain the observations.
CO isotologue emission lines show evidence for a global collapse of the hub with a corresponding large mass infall rate (i.e. ~10⁻³ M_⊙ yr⁻¹).

In conclusion, G345.88-1.10 is an intriguing source where the combined presence of an infrared dark fragment (H1) and associated mid-infrared bright bipolar cavities remains unexplained. The uniqueness of G345.88-1.10 strongly suggests that it traces a short and currently unknown phase of the high-mass star formation process.

The two main limitations of the study presented in this paper are the angular resolution of the available data (~20″–30″, i.e. 0.2 to 0.3 pc at the distance of G345.88-1.10) and the lack of shock tracers. ALMA observations at (sub)arcsecond resolution of the molecular outflow, fragments, cavity walls, and hub filaments would resolve a number of unsolved issues and place G345.88-1.10 in the larger context of high-mass star formation. Such observations would allow one to: i. determine the nature of the molecular outflow, i.e. whether it is made of a single outflow or multiple ones; ii. measure the mass of the core(s) that are driving the outflow(s) and e.g. place them in a luminosity/mass evolutionary diagram (Duarte-Cabral et al. 2013); iii. map the spatial distribution of shocked gas via the detection of Si0 emission; iv. map the velocity field of the hub at higher resolution, showing in more detail how the large-scale inflow of matter is connected to the smaller-scale core/disc.

Additionally, S0FIA observations of [C_II], [O_I], H₂, and high-J CO transitions would allow one to specifically look at the physics at play within the cavities, directly tackling the question of the heating mechanism. Lastly, sensitive arcsecond resolution radio continuum and RRL observations of the cavities with ATCA or a SKA precursor could allow the detection of an eventual radio jet at the origin of the cavities.

Acknowledgements

We thank the referee for insightful comments that significantly improved the content of this paper. We thank P. Clark and S. Mairs for helpful discussions on the cooling time and dust continuum variability of dense cores. This publication is based on data acquired with the Atacama Pathfinder Experiment (APEX) under program ID 097.F-9301(A). L.B. acknowledges support from the Erasmus+ program for his stay at Cardiff University. L.B. and N.S. acknowledge support by the French ANR and the German DFG through the GENESIS project (ANR-16-CE92-0035-01/DFG1591/2-1). N.P. acknowledges the support from STFC under the consolidated grant number ST/R005028/1. A.S. acknowledges the support of the Max Planck Society. This research made use of Astropy (https://www.astropy.org/; Astropy Collaboration 2013, 2018), a community-developed core Python package for Astronomy. This research made use of APLpy, an open-source plotting package for Python (Robitaille & Bressert 2012).

Appendix A Flux Determination

In order to calculate the nebulosity and core fluxes, we performed aperture photometry. For the nebulosity, we used two apertures of 25″ radius centred on each of the nebulosity lobes, encompassing most of the 70 μm emission. To correct for background emission, annuli were constructed around these apertures (see Fig. A.1). For the H1 fragment flux estimates, different aperture sizes were used depending on the wavelength and whether the core was unresolved. The aperture radius used for 160 μm and 250 μm, which resolve the fragment, was taken to be 13″ and centred on the H1 fragment as this corresponds to the fragment size. For the 350 μm and 500 μm data, the area was taken that corresponds to a single beam size (i.e. $\frac{π {FWHM}^{2}}{4 ln (2)}$ ${{\pi {\rm{FWH}}{{\rm{M}}^2}} \over {4{\rm{ln}}\left( 2 \right)}}$ ) centred on the H1 fragment. As for the nebulosity, annuli were constructed around the different apertures to correct for background emission, the annuli were taken to cover emission between 2 and 3 times the radius of the beamsize.

Fig. A.1

Apertures used to calculate the flux of the nebulosities in G345.88-1.10 (black). The annuli to correct for background emission are indicated with the dashed red circles.

Appendix B ¹²CO(2–1) Channel Maps

Fig. B.1 presents the channel maps of the ¹²CO(2–1) emission line. Focusing on the velocity ranges associated with the blueshifted wing (v_LSR < −25 km s⁻¹) and the redshifted wing (v_LSR > −15 km s⁻¹), it appears that there is weak high-velocity C0 emission over the full cavities and no evident velocity substructure in the outflow. This suggests that the observed protostellar outflow might be dominated by a single powerful outflow. Because of the limited resolution and sensitivity it is however challenging to be conclusive about this emission. Sensitive high-resolution observations will be able to shed more light on this.

Fig. B.1

¹²CO(2–1) channel maps over the full velocity range of the spectra with white Herschel column density contours highlighting the G345.88-1.10 hub. The colourbar was chosen to highlight the emission associated with the ¹²CO(2–1) wings, which explain why the figure is saturated in the figure between −25 km s⁻¹ and −15 km s⁻¹.

Appendix C The C¹⁸O(2–1) and ¹³CO(2–1) Emission

Fig. C.1 presents the velocity integrated C¹⁸O emission between −30 km s⁻¹ and −10 km s⁻¹ over the G345.88-1.10 hub. It is found that the C¹⁸O(2–1) emission broadly follows the morphology of the hub. However, the dust column density peaks of the fragments do not coincide perfectly with the C¹⁸O(2–1) peak emission. This could be due to several factors. First, CO could be depleted in the densest fragments. Second, strong temperature gradients could contribute to the line strength variations on its own. Further C¹⁸O observations of different transitions, which we currently lack, would help to constrain this question.

Some C¹⁸O(2–1) and ¹³CO(2–1) spectra from the hub are presented in Fig. C.1. In particular towards the cavities (spectra 5 & 6 in Fig. C.1), both the C¹⁸O(2–1) and ¹³CO(2–1) spectra show skewed line profiles. This presence of skewed line profiles in the cavities was also highlighted in Fig. 13. Towards the fragments, we observe blueshifted self-absorbed ¹³CO spectra towards fragment H1 and H2, a well-established signature of collapse signature which is widespread across the hub (see Sec. 8.2).

To obtain a centroid velocity and velocity dispersion maps, a single gaussian profile was fitted to the C¹⁸O(2–1) spectra. The resulting velocity field is presented in Fig. 14, showing organised kinematics with a velocity gradient peak centred on the H1 and H2 fragments. This strong velocity gradient towards the presumed centre of collapse is expected in the context of gravitational acceleration. In Fig. C.1, the velocity dispersion map is presented. It shows large spectral linewidths to the west of the hub, whereas the smallest linewidths are found to the east and north-east where the 10 pc long filament connects to the hub. Lastly, it can be noted that there is no observed increase in velocity dispersion towards the identified massive fragments which suggests they do not experience increased turbulent support compared to the rest of the hub.

Fig. C.1

(left) C¹⁸O(2–1) integrated intensity map from -30 km s⁻¹ to −10 km s⁻¹. Overplotted in black are the Herschel column density contours starting at $N_{H_{2}} = 2 \times 10^{22}$ ${{\rm{N}}_{{H_2}}} = 2 \times {10^{22}}$ cm⁻² with increments of 2×10²² cm⁻². The black crosses indicate the 4 extracted fragments from the column density map and the black rectangle indicates the region used to extract the spectra presented in Fig. 2. (right) The velocity dispersion map from fitting the Gaussian to the C¹⁸O(2–1) spectra, with overplotted in black the Herschel column density contours. The white crosses indicate the locations of the spectra presented below. (bottom) The ¹³CO(2–1) and C¹⁸O(2–1) spectra at the locations in the hub indicated in the velocity dispersion map. The spectra at the location of the cavities (spectra 5 & 6) show skewed spectra in both ¹³CO(2–1) and C¹⁸O(2–1).

Appendix D Outflow Momentum Rate Calculation

To determine the outflow momentum rate, we make use of the conservation of momentum rate along the outflow axis. For this exercise, an annulus with size ∆r allows to focus on a specific location along the outflow axis. The construction of the annulus thus allows to define a region along the outflow axis over which we work with conservation of momentum rate. Here, we chose ∆r = 28″ (~ 0.32 pc at the distance of G345.88-1.10), which is the APEX beam FWHM in ¹²CO(2–1). In a second step, a more complex shape was cut out from this annulus that is encompassed by the outflow contours, see Fig. D.1, to reduce noise from non-outflow material that is found in the annulus (Duarte-Cabral et al. 2013). From Fig. D.1, it is observed that the peak intensity is very close to the central source. The inner radius of the annulus around the outflow peak is thus unresolved. As a result, the annulus becomes effectively a circle with radius ∆r. Higher angular resolution observations will allow us to resolve the outflow structure, and provide a more accurate determination of its corresponding momentum rate.

The two annulus/circle for the outflow wings was constructed around the maximum of the outflow intensities to maximise the outflow rate since this should best trace the impact on the ambient cloud (Duarte-Cabral et al. 2013). Within the obtained areas, the momentum rate can be estimated for the outflows using: $F_{C O} \propto \int \frac{T (v) {(v - v_{0})}^{2} dv}{Δ r}$ ${F_{CO}} \propto \int {{{T\left( {\rm{v}} \right){{\left( {{\rm{v}} - {v_0}} \right)}^2}{\rm{dv}}} \over {\Delta {\rm{r}}}}}$ (D.1)

In this equation, T is the main beam temperature in Kelvin and v0 the velocity of the powering object which is assumed to be at −21 km s⁻¹. Further, a CO to H₂ conversion factor of 10⁴, a temperature of 20 K and a molecular weight of 2.33 g mole⁻¹ were used.

We include a correction factor for the opacity of the CO wings, given by $\frac{τ_{C O}}{1 - e^{- τ_{C O}}} = 3.5$ ${{{\tau _{CO}}} \over {1 - {e^{ - {\tau _{CO}}}}}} = 3.5$ (Cabrit & Bertout 1992). Outflow momentum rates for both the red- and blueshifted sides of the two outflows are given in Tab. 1.

Lastly, we also include a correction for the inclination of the outflow. The actual ∆r is given by ∆r/sin(i), and the real outflow velocity is given by v/cos(i). The correction for the inclination angle in the outflow momentum rate calculation is therefore given by $\frac{\sin (i)}{\cos^{2} (i)} = 32.7$ ${{\sin \left( {\rm{i}} \right)} \over {{{\cos }^2}\left( {\rm{i}} \right)}} = 32.7$ , which corresponds to the assumed inclination angle with the line of sight of 80°.

Fig. D.1

The circle (in black) with a FWHM radius encompasses the emission peak of the blue- and redshifted outflow. To reduce noise, the outer shape (fat line) of the outflow contours is then used to cut out the outflow inside the circle. This allows to focus only on the actual outflow located inside the circle and reduces the noise that is not part of the outflow. The white circle in the top left indicates the beamsize.

Appendix E Outflow Opening Angle

To determine the opening angle for the bipolar outflow cavities, a cone was constructed that connects the centre of the H1 fragment with the outer edge of the cavity. This is illustrated in Fig. E.1 for the 70 μm emission. This estimate of the outflow opening angle assumes axial symmetry which is uncertain seeing the proposed precession in the paper.

Fig. E.1

Illustration of the cones used to determine the opening angle of the 70 μm nebulosities.

Appendix F H30α non-detection

The spectra in Fig. F.1 present the H30α non-detection in both cavities. This provides 3σ upper limits of 63 mK and 110 mK for N1 and N2, respectively. Using the optically thin assumption for the H30α line, a size of 0.6 pc for the cavity and a typical electron temperature (T_e) of 7500 K for the Galactocentric distance of G345.88-1.10 (Quireza et al. 2006), it is then possible to calculate an upper limit for the average electron density. This can be done from the formula below (Gordon & Sorochenko 2002): $(\frac{T_{L}}{K}) ≃ 3076 {(\frac{T_{e}}{K})}^{- 3 / 2} (\frac{E M}{c m^{- 6} p c}) {(\frac{v_{0}}{G H z})}^{- 1} {(\frac{Δ V}{k m s^{- 1}})}^{- 1} Δ n M_{Δ n}$ $\left( {{{{T_L}} \over K}} \right) \simeq 3076{\left( {{{{T_e}} \over K}} \right)^{{{ - 3} \mathord{\left/ {\vphantom {{ - 3} 2}} \right. \kern-\nulldelimiterspace} 2}}}\left( {{{EM} \over {c{m^{ - 6}}pc}}} \right){\left( {{{{v_0}} \over {GHz}}} \right)^{ - 1}}{\left( {{{\Delta V} \over {km{s^{ - 1}}}}} \right)^{ - 1}}\Delta n{M_{\Delta n}}$ (F.1)

with T_L the line peak intensity, EM the emission measure, ν₀ the rest frequency of the line, ∆n = 1 and M_∆n = 0.19. The emission measure is approximated by $n_{e}^{2} L$ ${\rm{n}}_e^2{\rm{L}}$ , with n_e the electron density and L the size of the cavity in the line of sight. Plugging all values in this equation gives an electron density upper limit of n_e < 3.4×10² cm⁻³ and n_e < 4.6×10² cm⁻³ for the N1 and N2 cavities, respectively.

Fig. F.1

Average H30α spectra, sampled at 2 km s⁻¹ towards the N1 (top) and N2 (bottom) cavity. The vertical dashed lines indicates the expected central velocity of the line which remains undetected in both cavities.

Appendix G The Combination of a 70 μm Bright Cavity and 70 μm Quiet Core in the Radiative Transfer Simulations

We combined the gathered information on both the emission properties of the nebulosity and that of the core to check whether the radiative transfer simulations with nebulosities also correspond to the best core runs. To do this we took all the radiative transfer calculations for which a nebulosity has formed, and plotted the corresponding 70 μm nebulosity fluxes (normalised to the observed nebulosity flux) as a function of their associated $χ_{core}^{2}$ $\chi _{{\rm{core}}}^2$ (see Fig. G.1). First, we find that only a small number of runs with nebulosity have a $χ_{core}^{2} < 3$ $\chi _{{\rm{core}}}^2 < 3$ . Second, those good core runs are also the worst at reproducing the 70 μm nebulosity flux. In fact, these two conditions (i.e. a quiet mid-infrared core and a bright mid-infrared cavity) are conflicting, with the best models for each sitting at both ends of the point distribution in Fig. G.1. This important discrepancy and the fact that it appears impossible to even create radiative transfer nebulosities with more than 15% of the observed 70 μm flux in G345.88-1.10 strongly suggests that radiative heating from an embedded object is not responsible for the bright cavities in G345.88-1.10. This also justifies our approach not to apply a model optimizer and explore the parameter space by hand as it will search for a non existing solution.

As already pointed out in Sect. 7.4.2, radiative transfer calculations with a central source with L ≥ 4000 L_⊙ are not compatible with the observations of H1. Investigating the relation between the 70 μm flux of the core and the density of the models (see Fig. G.2), it is also found that a high density is necessary to reproduce a weak 70 μm core emission. At 500 L0, the central density n_c of the Plummer profile, or the characteristic density n₀ of the passive flared disc, needs to be larger than $n_{H_{2}} = 5 \times 10^{6}$ ${{\rm{n}}_{{H_2}}} = 5 \times {10^6}$ cm⁻³ to sufficiently suppress the 70 μm emission towards the core. For an embedded object with a luminosity of 2000 L_⊙ and more, the required density increases to values $n_{H_{2}} \geq 10^{7}$ ${{\rm{n}}_{{H_2}}} \ge {10^7}$ cm⁻³. However, as pointed out in Sec. 7.4.2, at 2000 L_⊙ there is still the problem that emission at longer wavelengths is not in agreement with what is observed towards the H1 fragment.

Fig. G.1

$χ_{c o r e}^{2}$ $\chi _{core}^2$ for the radiative transfer models that have nebulosities as a function of the 70 μm flux for the simulated nebulosities (normalised to the observed nebulosities). This shows that none of the radiative transfer simulations with 70 μm nebulosities fit with the core at longer wavelengths.

Fig. G.2

The 70 μm flux predicted by the radiative transfer calculations as a function of the central (n_c) and characteristic (n₀) density for all models with a luminosity of 500 L_⊙. The horizontal line indicates the upper limit for the 70 μm flux towards the H1 fragment.

Appendix H Calculating the Cooling Timescale for the H1 Fragment and the Cavities

The cooling time scale is given by $τ_{c o o l} = \frac{ρ u}{λ_{c o o l}}$ ${\tau _{cool}} = {{\rho {\rm{u}}} \over {{\lambda _{cool}}}}$ (H.1)

with ρ the gas density, u the specific internal energy and λ_cool the cooling rate per unit volume. The specific internal energy is given by $u = \frac{3}{2} \frac{k_{b} T_{g a s}}{μ m_{p}}$ ${\rm{u}} = {3 \over 2}{{{k_b}{{\rm{T}}_{gas}}} \over {\mu {{\rm{m}}_p}}}$ (H.2)

with k_b the Boltzmann constant, T_gas the kinetic gas temperature, μ the mean molecular mass (=2.33) and m_p the mass of a proton. The last term that has to be constrained is the cooling rate per unit volume (λ_coo1). If the gas is dense enough, typically $n_{H_{2}} > 10^{5}$ ${{\rm{n}}_{{H_2}}} > {10^5}$ cm⁻³ (Glover & Clark 2012), then it can be assumed that the cooling is driven by dust through gas-grain energy transfer. Given that the estimated average density in the H1 fragment is $n_{H_{2}} = 3.2 \times 10^{5}$ ${{\rm{n}}_{{H_2}}} = 3.2 \times {10^5}$ cm⁻³, we can use this condition (Glover & Clark 2012). When the cooling is dominated by gas-grain energy transfer the cooling rate per unit volume is given by $λ_{c o o l} = 4 π ρ D \int_{0}^{\infty} κ_{v} B_{v} (T) d v$ ${\lambda _{cool}} = 4\pi \rho {\rm{D}}\int_0^\infty {{\kappa _v}{{\rm{B}}_v}\left( {\rm{T}} \right){\rm{d}}v}$ (H.3)

where κ_ν is the dust opacity, D the dust to gas mass fraction and B_ν the Planck function. Using Ossenkopf & Henning (1994) this results in $λ_{c o o l} = 4.68 \times 10^{- 31} {nT}_{d u s t}^{6} (in erg s^{- 1} {cm}^{- 3})$ ${\lambda _{cool}} = 4.68 \times {10^{ - 31}}{\rm{nT}}_{dust}^6\left( {{\rm{in}}\,{\rm{erg}}\,{{\rm{s}}^{ - 1}}\,{\rm{c}}{{\rm{m}}^{ - 3}}} \right)$ (H.4)

with T_dust being the dust temperature and n being the number density of the gas for which we assume it is fully molecular in the H1 fragment.

Table H.1

Cooling timescales for the H1 fragment for a grid of dust and gas temperatures.

We considered several input parameters for the heated gas in the fragment ranging from 20 K to 100 K and for T_dust we considered that the heated dust could reach 20 K, 30 K and 40 K. The resulting cooling timescales are given in Tab. H.1 and are typically of the order of a few days to a few years for these conditions. Assuming even higher dust temperature values would only further decrease the cooling timescale.

For the cavities we get typical densities $n_{H_{2}} = 1.0 - 2.5 \times 10^{4}$ ${{\rm{n}}_{{H_2}}} = 1.0 - 2.5 \times {10^4}$ cm⁻³ based on the Herschel column density and cavity size. From Goldsmith & Langer (1978) we find that for a typical density $n_{H_{2}} = 10^{4}$ ${{\rm{n}}_{{H_2}}} = {10^4}$ cm⁻³, λ_coo1 is given by $λ_{c o o l} = 1.5 \times 10^{- 26} T^{2.7} (in erg s^{- 1} {cm}^{- 3})$ ${\lambda _{cool}} = 1.5 \times {10^{ - 26}}{{\rm{T}}^{2.7}}\left( {{\rm{in}}\,{\rm{erg}}\,{{\rm{s}}^{ - 1}}\,{\rm{c}}{{\rm{m}}^{ - 3}}} \right)$ (H.5)

where T is the kinetic temperature in the region. Using the same average kinetic temperatures of 20–100 K in the cavities, we obtain λ_coo1 between 4.9×10⁻²³ − 3.8×10⁻²¹ erg cm⁻³ s⁻¹ whichis consistent with cooling rates found in the simulations by Glover & Clark (2012) for these densities. Using these values, we obtain estimated cooling timescales between 1.7×10³ and 2.7×10⁴ yr for the cavities. Note that the lower end of these estimated cooling timescales is of the same order of magnitude as the propagation timescale for the jet through the cavity. This could thus be consistent with the observation that the full cavity is bright and not a specific subregion, Because of the lower densities, the cavities thus have significantly longer cooling timescales than the dense fragment.

References

Adams, J.D., Herter, T.L., Hora, J.L., et al. 2015, ApJ, 814, 54 [NASA ADS] [CrossRef] [Google Scholar]
Ahmadi, A., Beuther, H., Mottram, J.C., et al. 2018, A&A, 618, A46 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Anderson, M., Peretto, N., Ragan, S.E., et al. 2021, MNRAS, 508, 2964 [NASA ADS] [CrossRef] [Google Scholar]
André, P., Ward-Thompson, D., & Barsony, M. 2000, Protostars and Planets IV (Tucson: University of Arizona Press), 59 [Google Scholar]
André, P., Di Francesco, J., Ward-Thompson, D., et al. 2014, in Protostars and Planets VI, eds. H. Beuther, R.S. Klessen, C.P. Dullemond, & T. Henning (Tucson: University of Arizona Press), 27 [Google Scholar]
Anglada, G., Rodriguez, L.F., & Carrasco-Gonzâlez, C. 2018, A&ARv, 26, 3 [NASA ADS] [CrossRef] [Google Scholar]
Arce, H.G., & Sargent, A.I. 2006, ApJ, 646, 1070 [NASA ADS] [CrossRef] [Google Scholar]
Arce, H.G., Shepherd, D., Gueth, F., et al. 2007, Protostars and Planets V (Tucson: University of Arizona Press), 245 [Google Scholar]
Arce, H.G., Borkin, M.A., Goodman, A.A., Pineda, J.E., & Beaumont, C.N. 2011, ApJ, 742, 105 [NASA ADS] [CrossRef] [Google Scholar]
Arzoumanian, D., Russeil, D., Zavagno, A., et al. 2022, A&A, 660, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Astropy Collaboration (Robitaille, T.P., et al.) 2013, A&A, 558, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Astropy Collaboration (Price-Whelan, A.M., et al.) 2018, AJ, 156, 123 [NASA ADS] [CrossRef] [Google Scholar]
Avison, A., Fuller, G.A., Peretto, N., et al. 2021, A&A, 645, A142 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Bachiller, R. 1996, ARA&A, 34, 111 [NASA ADS] [CrossRef] [Google Scholar]
Bally, J. 2016, ARA&A, 54, 491 [Google Scholar]
Battersby, C., Bally, J., Ginsburg, A., et al. 2011, A&A, 535, A128 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Beltrán, M.T., & de Wit, W.J. 2016, A&ARv, 24, 6 [CrossRef] [Google Scholar]
Bendo, G.J., Griffin, M.J., Bock, J.J., et al. 2013, MNRAS, 433, 3062 [NASA ADS] [CrossRef] [Google Scholar]
Benjamin, R.A., Churchwell, E., Babler, B.L., et al. 2003, PASP, 115, 953 [NASA ADS] [CrossRef] [Google Scholar]
Berger, M.J., & Colella, P. 1989, J. Comput. Phys., 82, 64 [CrossRef] [Google Scholar]
Berger, M.J., & Oliger, J. 1984, J. Comput. Phys., 53, 484 [NASA ADS] [CrossRef] [Google Scholar]
Beuther, H., Henning, T., Linz, H., et al. 2015a, A&A, 581, A119 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Beuther, H., Ragan, S.E., Johnston, K., et al. 2015b, A&A, 584, A67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Beuther, H., Mottram, J.C., Ahmadi, A., et al. 2018, A&A, 617, A100 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Bjorkman, J.E., & Wood, K. 2001, ApJ, 554, 615 [NASA ADS] [CrossRef] [Google Scholar]
Bock, D.C.-J., Large, M.I., & Sadler, E.M. 1999, AJ, 117, 1578 [NASA ADS] [CrossRef] [Google Scholar]
Bonne, L., Bontemps, S., Schneider, N., et al. 2020, A&A, 644, A27 [EDP Sciences] [Google Scholar]
Bonne, L., Schneider, N., Garcia, P., et al. 2022, ApJ, accepted [arXiv:2287.86479] [Google Scholar]
Bontemps, S., Andre, P., Terebey, S., & Cabrit, S. 1996, A&A, 311, 858 [NASA ADS] [Google Scholar]
Bontemps, S., Motte, F., Csengeri, T., & Schneider, N. 2010, A&A, 524, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Borissova, J., Bonatto, C., Kurtev, R., et al. 2011, A&A, 532, A131 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Cabrit, S., & Bertout, C. 1986, ApJ, 307, 313 [NASA ADS] [CrossRef] [Google Scholar]
Cabrit, S., & Bertout, C. 1992, A&A, 261, 274 [Google Scholar]
Carey, S.J., Noriega-Crespo, A., Mizuno, D.R., et al. 2009, PASP, 121, 76 [NASA ADS] [CrossRef] [Google Scholar]
Carrasco-Gonzâlez, C., Rodriguez, L.F., Anglada, G., et al. 2010, Science, 330, 1209 [CrossRef] [Google Scholar]
Carrasco-Gonzâlez, C., Torrelles, J.M., Cantö, J., et al. 2015, Science, 348, 114 [CrossRef] [Google Scholar]
Cesaroni, R., Sânchez-Monge, Ä., Beltrân, M.T., et al. 2017, A&A, 602, A59 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Cesaroni, R., Moscadelli, L., Neri, R., et al. 2018, A&A, 612, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Comerön, F., Schneider, N., Djupvik, A.A., & Schnugg, C. 2018, A&A, 615, A2 [CrossRef] [EDP Sciences] [Google Scholar]
Condon, J.J., Griffith, M.R., & Wright, A.E. 1993, AJ, 106, 1095 [NASA ADS] [CrossRef] [Google Scholar]
Csengeri, T., Bontemps, S., Wyrowski, F., et al. 2017a, A&A, 601, A60 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Csengeri, T., Bontemps, S., Wyrowski, F., et al. 2017b, A&A, 600, L10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Csengeri, T., Bontemps, S., Wyrowski, F., et al. 2018, A&A, 617, A89 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Deharveng, L., Zavagno, A., Samal, M.R., et al. 2015, A&A, 582, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
de Valon, A., Dougados, C., Cabrit, S., et al. 2020, A&A, 634, L12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
De Vries, C.H., & Myers, P.C. 2005, ApJ, 620, 800 [NASA ADS] [CrossRef] [Google Scholar]
Dhabal, A., Mundy, L.G., Rizzo, M.J., Storm, S., & Teuben, P. 2018, ApJ, 853, 169 [NASA ADS] [CrossRef] [Google Scholar]
di Francesco, J., Evans, N.J.I., Caselli, P., et al. 2007, in Protostars and Planets V, eds. B. Reipurth, D. Jewitt, & K. Keil (Tucson: University of Arizona Press), 17 [Google Scholar]
Draine, B.T. 2011, Physics of the Interstellar and Intergalactic Medium (Priceton: Princeton University Press) [Google Scholar]
Draine, B.T., & Li, A. 2007, ApJ, 657, 810 [NASA ADS] [CrossRef] [Google Scholar]
Duarte-Cabral, A., Chrysostomou, A., Peretto, N., et al. 2012, A&A, 543, A140 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Duarte-Cabral, A., Bontemps, S., Motte, F., et al. 2013, A&A, 558, A125 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Duarte-Cabral, A., Bontemps, S., Motte, F., et al. 2014, A&A, 570, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Dullemond, C.P. 2012, Astrophysics Source Code Library [record ascl:1282.815] [Google Scholar]
Ellerbroek, L.E., Bik, A., Kaper, L., et al. 2013, A&A, 558, A102 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Fazio, G.G., Hora, J.L., Allen, L.E., et al. 2004, ApJS, 154, 10 [NASA ADS] [CrossRef] [Google Scholar]
Frank, A., Ray, T.P., Cabrit, S., et al. 2014, in Protostars and Planets VI, eds. H. Beuther, R.S. Klessen, C.P. Dullemond, & T. Henning (Tucson: University of Arizona Press), 451 [Google Scholar]
Froebrich, D., Scholz, A., & Raftery, C.L. 2007, MNRAS, 374, 399 [NASA ADS] [CrossRef] [Google Scholar]
Glover, S.C.O., & Clark, P.C. 2012, MNRAS, 421, 9 [Google Scholar]
Goldsmith, P.F., & Langer, W.D. 1978, ApJ, 222, 881 [NASA ADS] [CrossRef] [Google Scholar]
Gömez, G.C., & Vâzquez-Semadeni, E. 2014, ApJ, 791, 124 [CrossRef] [Google Scholar]
Gordon, M.A., & Sorochenko, R.L. 2002, Radio Recombination Lines. Their Physics and Astronomical Applications (Berlin:Springer), 282 [Google Scholar]
Griffin, M.J., Abergel, A., Abreu, A., et al. 2010, A&A, 518, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Griffin, M.J., North, C.E., Schulz, B., et al. 2013, MNRAS, 434, 992 [NASA ADS] [CrossRef] [Google Scholar]
Griffith, M.R., & Wright, A.E. 1993, AJ, 105, 1666 [NASA ADS] [CrossRef] [Google Scholar]
Gueth, F., Guilloteau, S., & Bachiller, R. 1998, A&A, 333, 287 [NASA ADS] [Google Scholar]
Gusdorf, A., Anderl, S., Lefloch, B., et al. 2017, A&A, 602, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Güsten, R., Nyman, L.Â., Schilke, P., et al. 2006, A&A, 454, L13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Habel, N.M., Megeath, S.T., Booker, J.J., et al. 2021, ApJ, 911, 153 [NASA ADS] [CrossRef] [Google Scholar]
Hacar, A., Alves, J., Tafalla, M., & Goicoechea, J.R. 2017, A&A, 602, L2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Hartmann, L., & Burkert, A. 2007, ApJ, 654, 988 [NASA ADS] [CrossRef] [Google Scholar]
Hartmann, L., Herczeg, G., & Calvet, N. 2016, ARA&A, 54, 135 [Google Scholar]
Henshaw, J.D., Caselli, P., Fontani, F., et al. 2013, MNRAS, 428, 3425 [NASA ADS] [CrossRef] [Google Scholar]
Henshaw, J.D., Caselli, P., Fontani, F., Jiménez-Serra, I., & Tan, J.C. 2014, MNRAS, 440, 2860 [CrossRef] [Google Scholar]
Hildebrand, R.H. 1983, QJRAS, 24, 267 [Google Scholar]
Hill, T., Motte, F., Didelon, P., et al. 2011, A&A, 533, A94 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Hunter, T.R., Brogan, C.L., De Buizer, J.M., et al. 2021, ApJ, 912, L17 [CrossRef] [Google Scholar]
Ilee, J.D., Cyganowski, C.J., Nazari, P., et al. 2016, MNRAS, 462, 4386 [NASA ADS] [CrossRef] [Google Scholar]
Ilee, J.D., Cyganowski, C.J., Brogan, C.L., et al. 2018, ApJ, 869, L24 [NASA ADS] [CrossRef] [Google Scholar]
Inoue, T., Hennebelle, P., Fukui, Y., et al. 2018, PASJ, 70, S53 [NASA ADS] [Google Scholar]
Jackson, J.M., Whitaker, J.S., Rathborne, J.M., et al. 2019, ApJ, 870, 5 [Google Scholar]
Johnston, K.G., Keto, E., Robitaille, T.P., & Wood, K. 2011, MNRAS, 415, 2953 [NASA ADS] [CrossRef] [Google Scholar]
Johnston, K.G., Robitaille, T.P., Beuther, H., et al. 2015, ApJ, 813, L19 [NASA ADS] [CrossRef] [Google Scholar]
Johnstone, D., Hendricks, B., Herczeg, G.J., & Bruderer, S. 2013, ApJ, 765, 133 [NASA ADS] [CrossRef] [Google Scholar]
Kalcheva, I.E., Hoare, M.G., Urquhart, J.S., et al. 2018, A&A, 615, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Kirk, H., Myers, P.C., Bourke, T.L., et al. 2013, ApJ, 766, 115 [NASA ADS] [CrossRef] [Google Scholar]
Krumholz, M.R., & Tan, J.C. 2007, ApJ, 654, 304 [NASA ADS] [CrossRef] [Google Scholar]
Kuiper, R., Klahr, H., Beuther, H., & Henning, T. 2010, ApJ, 722, 1556 [NASA ADS] [CrossRef] [Google Scholar]
Kuiper, R., Turner, N.J., & Yorke, H.W. 2016, ApJ, 832, 40 [NASA ADS] [CrossRef] [Google Scholar]
Larson, R.B., & Starrfield, S. 1971, A&A, 13, 190 [NASA ADS] [Google Scholar]
Lee, Y.-H., Johnstone, D., Lee, J.-E., et al. 2021, ApJ, 920, 119 [NASA ADS] [CrossRef] [Google Scholar]
Lehmann, A., Godard, B., Pineau des Forêts, G., & Falgarone, E. 2020, A&A, 643, A101 [EDP Sciences] [Google Scholar]
Lehmann, A., Godard, B., Pineau des Forêts, G., Vidal-Garcia, A., & Falgarone, E. 2022, A&A, 658, A165 [CrossRef] [EDP Sciences] [Google Scholar]
Li, Z.Y., Banerjee, R., Pudritz, R.E., et al. 2014, in Protostars and Planets VI, eds. H. Beuther, R.S. Klessen, C.P. Dullemond, & T. Henning (Tucson: University of Arizona Press), 173 [Google Scholar]
Louvet, F., Dougados, C., Cabrit, S., et al. 2018, A&A, 618, A120 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Louvet, F., Neupane, S., Garay, G., et al. 2019, A&A, 622, A99 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Lumsden, S.L., Hoare, M.G., Urquhart, J.S., et al. 2013, ApJS, 208, 11 [NASA ADS] [CrossRef] [Google Scholar]
Mairs, S., Johnstone, D., Kirk, H., et al. 2017, ApJ, 849, 107 [NASA ADS] [CrossRef] [Google Scholar]
Mardones, D., Myers, P.C., Tafalla, M., et al. 1997, ApJ, 489, 719 [NASA ADS] [CrossRef] [Google Scholar]
Mathis, J.S., Mezger, P.G., & Panagia, N. 1983, A&A, 128, 212 [Google Scholar]
Mauch, T., Murphy, T., Buttery, H.J., et al. 2003, MNRAS, 342, 1117 [NASA ADS] [CrossRef] [Google Scholar]
Minier, V., Tremblin, P., Hill, T., et al. 2013, A&A, 550, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Molinari, S., Swinyard, B., Bally, J., et al. 2010, A&A, 518, L100 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Molinari, S., Schisano, E., Elia, D., et al. 2016, A&A, 591, A149 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Motte, F., Bontemps, S., Schilke, P., et al. 2007, A&A, 476, 1243 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Motte, F., Zavagno, A., Bontemps, S., et al. 2010, A&A, 518, L77 [CrossRef] [EDP Sciences] [Google Scholar]
Motte, F., Bontemps, S., & Louvet, F. 2018, ARA&A, 56, 41 [Google Scholar]
Myers, P.C., Mardones, D., Tafalla, M., Williams, J.P., & Wilner, D.J. 1996, ApJ, 465, L133 [NASA ADS] [CrossRef] [Google Scholar]
Orkisz, J.H., Peretto, N., Pety, J., et al. 2019, A&A, 624, A113 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Ossenkopf, V., & Henning, T. 1994, A&A, 291, 943 [NASA ADS] [Google Scholar]
Palau, A., Fuente, A., Girart, J.M., et al. 2013, ApJ, 762, 120 [NASA ADS] [CrossRef] [Google Scholar]
Peretto, N., André, P., & Belloche, A. 2006, A&A, 445, 979 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Peretto, N., Hennebelle, P., & André, P. 2007, A&A, 464, 983 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Peretto, N., Fuller, G.A., Plume, R., et al. 2010, A&A, 518, L98 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Peretto, N., Fuller, G.A., Duarte-Cabral, A., et al. 2013, A&A, 555, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Peretto, N., Fuller, G.A., André, P., et al. 2014, A&A, 561, A83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Peretto, N., Lenfestey, C., Fuller, G.A., et al. 2016, A&A, 590, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Peretto, N., Rigby, A., André, P., et al. 2020, MNRAS, 496, 3482 [NASA ADS] [CrossRef] [Google Scholar]
Pineda, J.E., Arce, H.G., Schnee, S., et al. 2011, ApJ, 743, 201 [NASA ADS] [CrossRef] [Google Scholar]
Pineda, J.L., Langer, W.D., Velusamy, T., & Goldsmith, P.F. 2013, A&A, 554, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Poglitsch, A., Waelkens, C., Geis, N., et al. 2010, A&A, 518, L2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Qin, S.-L., Schilke, P., Rolffs, R., et al. 2011, A&A, 530, L9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Quireza, C., Rood, R.T., Bania, T.M., Balser, D.S., & Maciel, W.J. 2006, ApJ, 653, 1226 [NASA ADS] [CrossRef] [Google Scholar]
Reid, M.J., Menten, K.M., Brunthaler, A., et al. 2014, ApJ, 783, 130 [NASA ADS] [CrossRef] [Google Scholar]
Rigby, A.J., Peretto, N., Adam, R., et al. 2021, MNRAS, 502, 4576 [NASA ADS] [CrossRef] [Google Scholar]
Robitaille, T., & Bressert, E. 2012, Astrophysics Source Code Library [record ascl:1288.817] [Google Scholar]
Rodriguez-Kamenetzky, A., Carrasco-Gonzâlez, C., Araudo, A., et al. 2017, ApJ, 851, 16 [NASA ADS] [CrossRef] [Google Scholar]
Rosen, A.L., & Krumholz, M.R. 2020, AJ, 160, 78 [NASA ADS] [CrossRef] [Google Scholar]
Rosolowsky, E.W., Pineda, J.E., Kauffmann, J., & Goodman, A.A. 2008, ApJ, 679, 1338 [CrossRef] [Google Scholar]
Russeil, D., Zavagno, A., Motte, F., et al. 2010, A&A, 515, A55 [CrossRef] [EDP Sciences] [Google Scholar]
Sadavoy, S.I., Di Francesco, J., Johnstone, D., et al. 2013, ApJ, 767, 126 [NASA ADS] [CrossRef] [Google Scholar]
Samal, M.R., Deharveng, L., Zavagno, A., et al. 2018, A&A, 617, A67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Sandell, G., Wright, M., Güsten, R., et al. 2020, ApJ, 904, 139 [NASA ADS] [CrossRef] [Google Scholar]
Sanhueza, P., Contreras, Y., Wu, B., et al. 2019, ApJ, 886, 102 [Google Scholar]
Sault, R.J., Teuben, P.J., & Wright, M.C.H. 1995, ASP Conf. Ser., 77, 433 [NASA ADS] [Google Scholar]
Schisano, E., Rygl, K.L.J., Molinari, S., et al. 2014, ApJ, 791, 27 [NASA ADS] [CrossRef] [Google Scholar]
Schmiedeke, A., Schilke, P., Möller, T., et al. 2016, A&A, 588, A143 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Schneider, N., Csengeri, T., Bontemps, S., et al. 2010, A&A, 520, A49 [CrossRef] [EDP Sciences] [Google Scholar]
Schneider, N., Csengeri, T., Hennemann, M., et al. 2012, A&A, 540, L11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Schneider, N., Bontemps, S., Girichidis, P., et al. 2015a, MNRAS, 453, L41 [CrossRef] [Google Scholar]
Schneider, N., Csengeri, T., Klessen, R.S., et al. 2015b, A&A, 578, A29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Schneider, N., Röllig, M., Simon, R., et al. 2018, A&A, 617, A45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Schuller, F., Csengeri, T., Urquhart, J.S., et al. 2017, A&A, 601, A124 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Schuller, F., Urquhart, J.S., Csengeri, T., et al. 2021, MNRAS, 500, 3064 [Google Scholar]
Skrutskie, M.F., Cutri, R.M., Stiening, R., et al. 2006, AJ, 131, 1163 [NASA ADS] [CrossRef] [Google Scholar]
Smith, R.J., Shetty, R., Beuther, H., Klessen, R.S., & Bonnell, I.A. 2013, ApJ, 771, 24 [NASA ADS] [CrossRef] [Google Scholar]
Stecklum, B., Wolf, V., Linz, H., et al. 2021, A&A, 646, A161 [EDP Sciences] [Google Scholar]
Tabone, B., Cabrit, S., Bianchi, E., et al. 2017, A&A, 607, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Tafalla, M. 2017, Mem. Soc. Astron. It., 88, 619 [NASA ADS] [Google Scholar]
Tafalla, M., Myers, P.C., Caselli, P., & Walmsley, C.M. 2004, A&A, 416, 191 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Tan, J.C., Kong, S., Zhang, Y., et al. 2016, ApJ, 821, L3 [NASA ADS] [CrossRef] [Google Scholar]
Tigé, J., Motte, F., Russeil, D., et al. 2017, A&A, 602, A77 [Google Scholar]
Torrelles, J.M., Patel, N.A., Curiel, S., et al. 2011, MNRAS, 410, 627 [NASA ADS] [CrossRef] [Google Scholar]
Traficante, A., Calzoletti, L., Veneziani, M., et al. 2011, MNRAS, 416, 2932 [Google Scholar]
Vassilev, V., Meledin, D., Lapkin, I., et al. 2008, A&A, 490, 1157 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Ward-Thompson, D., André, P., Crutcher, R., et al. 2007, in Protostars and Planets V, eds. B. Reipurth, D. Jewitt, & K. Keil (Tucson: University of Arizona Press), 33 [Google Scholar]
Watkins, E.J., Peretto, N., Marsh, K., & Fuller, G.A. 2019, A&A, 628, A21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Wenger, T.V., Balser, D.S., Anderson, L.D., & Bania, T.M. 2018, ApJ, 856, 52 [NASA ADS] [CrossRef] [Google Scholar]
Williams, G.M., Peretto, N., Avison, A., Duarte-Cabral, A., & Fuller, G.A. 2018, A&A, 613, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Wolfire, M.G., & Cassinelli, J.P. 1987, ApJ, 319, 850 [NASA ADS] [CrossRef] [Google Scholar]
Wright, E.L., Eisenhardt, P.R.M., Mainzer, A.K., et al. 2010, AJ, 140, 1868 [NASA ADS] [CrossRef] [Google Scholar]
Yorke, H.W., & Bodenheimer, P. 1999, ApJ, 525, 330 [NASA ADS] [CrossRef] [Google Scholar]
Yorke, H.W., & Sonnhalter, C. 2002, ApJ, 569, 846 [NASA ADS] [CrossRef] [Google Scholar]
Zhang, Y., & Tan, J.C. 2018, ApJ, 853, 18 [NASA ADS] [CrossRef] [Google Scholar]
Zinnecker, H., & Yorke, H.W. 2007, ARA&A, 45, 481 [NASA ADS] [CrossRef] [Google Scholar]

¹

http://vialactea.iaps.inaf.it/

²

https://www.iram.fr/IRAMFR/GILDAS/doc/html/class-html/class.html

³

http://www.treywenger.com/kd/index.php

⁴

https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.brentq.html

⁵

http://www.dendrograms.org/

All Tables

Table 1

Identified fragment properties.

	Fig. 2 ¹²CO(2–1), ¹³CO(2–1), and C¹⁸O(2–1) emission averaged over the extent of the hub.
In the text

	Fig. 5 Mid-infrared images from IRAC (3.6, 4.5, 5.8 and 8 μm), WISE (22 μm) and Herschel (160, 250, 350, 500 μm) of the G345.88-1.10 hub with the overlaid contour indicating the area of the 70 μm nebulosities (here defined by being brighter than 2000 MJy sr⁻¹).
In the text

	Fig. 6 Spectral energy distributions obtained from Herschel, WISE and IRAC images, for the two nebulosities N1 and N2, and the H1 fragment.
In the text

	Fig. 8 Radio continuum image of G345.88-1.10 at 843 MHz from the SUMSS survey (Bock et al. 1999). The 70 μm emission observed towards the bipolar cavity is highlighted with black contours (2000, 4000, 6000, 8000 and 10 000 MJy sr⁻¹).
In the text

	Fig. 11 $χ_{core}^{2}$ $\chi _{{\mathop{\rm core}\nolimits} }^2$ as a function of the central source luminosity for all models with $χ_{core}^{2} < 3$ $\chi _{{\mathop{\rm core}\nolimits} }^2 < 3$ . Model 1: red circle; Model 2: blue triangle; Model 3: green square.
In the text

	Fig. 12 70 μm flux distributions for the simulated cores with $χ_{core}^{2} \leq 2$ $\chi _{{\mathop{\rm core}\nolimits} }^2 \le 2$ for an internal luminosity of 500 L_⊙. The vertical dashed line indicates the upper limit of 5 Jy observed towards G345.88-1.10.
In the text

	Fig. 13 ¹²CO(2–1), ¹³CO(2–1) and C¹⁸O(2–1) spectra extracted in the N1 (top) and N2 (bottom) cavities. It is observed that the high-velocity outflow wing, seen in ¹²CO(2–1), are at the opposite side of the spectrum with respect to the lower-velocity tails that are also seen in ¹³CO(2–1) and C¹⁸O(2–1) on top of ¹²CO(2–1).
In the text

	Fig. A.1 Apertures used to calculate the flux of the nebulosities in G345.88-1.10 (black). The annuli to correct for background emission are indicated with the dashed red circles.
In the text

	Fig. B.1 ¹²CO(2–1) channel maps over the full velocity range of the spectra with white Herschel column density contours highlighting the G345.88-1.10 hub. The colourbar was chosen to highlight the emission associated with the ¹²CO(2–1) wings, which explain why the figure is saturated in the figure between −25 km s⁻¹ and −15 km s⁻¹.
In the text

	Fig. E.1 Illustration of the cones used to determine the opening angle of the 70 μm nebulosities.
In the text

	Fig. F.1 Average H30α spectra, sampled at 2 km s⁻¹ towards the N1 (top) and N2 (bottom) cavity. The vertical dashed lines indicates the expected central velocity of the line which remains undetected in both cavities.
In the text

	Fig. G.1 $χ_{c o r e}^{2}$ $\chi _{core}^2$ for the radiative transfer models that have nebulosities as a function of the 70 μm flux for the simulated nebulosities (normalised to the observed nebulosities). This shows that none of the radiative transfer simulations with 70 μm nebulosities fit with the core at longer wavelengths.
In the text

	Fig. G.2 The 70 μm flux predicted by the radiative transfer calculations as a function of the central (n_c) and characteristic (n₀) density for all models with a luminosity of 500 L_⊙. The horizontal line indicates the upper limit for the 70 μm flux towards the H1 fragment.
In the text