A new phase of massive star formation? A luminous outflow cavity centred on an infrared quiet core

We present APEX, infrared and radio continuum observations 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. At a distance of 2.26$^{+0.30}_{-0.21}$ 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 densest fragment (with M=210 M$_{\odot}$, R$_{\rm{eff}}=0.14$ pc) sits at the centre of a wide (opening angle of $\sim$ 90$\pm$15$^{o}$) bipolar nebulosity. $^{12}$CO(2-1) observations show that these infrared-bright nebulosities are spatially associated with a powerful molecular outflow from the central fragment. Negligible radio continuum and no H30$\alpha$ emission is detected towards the cavities, seemingly excluding that ionising radiation drives the evolution of the cavities. Furthermore, radiative transfer simulations are unable to reproduce the observed combination of a low-luminosity ($\lesssim$ 500 L$_{\odot}$) central source and a surrounding high-luminosity ($\sim 4000$ L$_{\odot}$) mid-infrared-bright bipolar cavity. This suggests that radiative heating from a central protostar cannot be responsible for the illumination of the outflow cavities. To our knowledge, this is the first reported object of this type. The rarity of objects like G345.88-1.10 is likely related to a very short phase in the massive star and/or cluster formation process that was so far unidentified. 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.


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 ). 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. 2007Russeil 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(Peretto et al. , 2013(Peretto et al. , 2014Schneider et al. 2010Schneider et al. , 2015aKirk 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 lowmass 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 socalled 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 highangular 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 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.

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(Schuller et al. , 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).

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 1 (Molinari et al. 2010(Molinari et al. , 2016Traficante 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.
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). 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.

APEX
In march 2016, the 12 CO(2-1), 13 CO(2-1) and C 18 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 12 CO(2-1) simultaneously also covers the H30α hydrogen recombination line at 231.90 GHz, which is a tracer of H II regions. This data has an angular resolution of ∼28 . The 13 CO(2-1) and C 18 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 2 with a channel width of 0.3 km s −1 for all lines. The noise, estimated as the line standard deviation within signal-free channels, is ∼0.2 K for 12 CO(2-1) and H30α. For 13 CO(2-1) and C 18 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 α (2000) = 17 h 11 m 43.9 s , δ (2000) = −41 • 18 25.5 . The 12 CO(2-1) and H30α map covers an area on the sky of 4 × 4 while the 13 CO(2-1) and C 18 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). Our APEX data show (see Fig. 2) that the velocity of G345.88-1.10 is centred at −21 km s −1 . This results in a heliocentric kinematic near distance of 2.26 +0.30 −0.21 kpc, where the uncertainties are given by the Monte Carlo calculations 3 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.24 −0.48 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 5 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 5 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 −1 and −6 km s −1 that we identify in Fig. 2 would be foreground clouds. These two clouds overlap in velocity with the observed CO 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  The H 2 column density contours are overplotted in white starting at N H 2 = 2 × 10 22 cm −2 with increments of 2 × 10 22 cm −2 . 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 Herschel dust temperature centred on the infrared nebulosity. All symbols and notations are identical to those presented in the left-hand-side panel.
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. Figure 1 presents the H 2 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 Herschel. First, we convolved the 160 µm image to the resolution of the 250 µm image using a Gaussian kernel of FWHM = θ 2 250 − θ 2 160 = 13.4 . Then, the temperature at every pixel was determined from the ratio of the 160 and 250 µm flux density using

H 2 column density and dust temperature mapping
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 4 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 2 column density image can then be obtained via the following equation: where κ ν 250 = 0.012 cm 2 g −1 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). 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.

Identifying filaments and fragments
In order to identify sub-structures within the column density map we used the astrodendro package 5 (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 cm −2 , to focus on the dense gas only, and is constructed using steps of N H 2 = 2.5 × 10 21 cm −2 . With these  (6): averaged H 2 column density. Errors on the size are estimated based on the distance uncertainty. For the column density an error of 50% is assumed to take into account uncertainties related to the dust opacity (Ossenkopf & Henning 1994), data calibration Bendo et al. 2013) and dust to gas conversion. The error on the mass includes both the 50% as well as the distance uncertainty.
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 +50 −30 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.

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 2 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, Hband 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.

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. Only 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.2 −1.1 × 10 3 L and 2.3 +1.3 −1.2 × 10 3 L , respectively. The luminosity of the H1 fragment was found to be 51 +29 −27 L only. The luminosity of the entire region (N1 + N2 + H1) is 4.4 × 10 3 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.

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.  <5.0 22 ± 4 27 ± 5 22 ± 4 15 ± 3 N1 0.25 ± 0.05 0.51 ± 0.10 0.38 ± 0.08 0.27 ± 0.05 3.1 ± 0.6 8.0 ± 1.6 4.0 ± 0.8 9.9 ± 2.0 (1.5 ± 0.3) × 10 2 53 ± 11 N2 0.21±0.04 0.40 ± 0.08 0.34 ± 0.07 0.31 ± 0.06 3.2 ± 0.6 8.2 ± 1.6 4.3 ± 0.9 9.0 ± 1.8 (1.8 ± 0.4) × 10 2 71 ± 14 Notes. An error of 20% is considered for all fluxes, to take into account uncertainty on the calibration and aperture. The presence of outflows can be verified by the inspection of the 12 CO(2-1) data. In fact, in Fig. 2, one can already notice the presence of high-velocity wings in the average 12 CO(2-1) spectrum of the hub, a clear signature of outflowing gas. By integrating these 12 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 12 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 18 O(2-1) centroid velocity (v c 18 o ) and velocity dispersion (σ c 18 o ) were first derived for each pixel (see Appendix C). After visual inspection, we defined the 12 CO outflow velocity ranges as Finally, as can be seen in Fig. 2, contaminating emission is present around v lsr = −5 km s −1 km s −1 and v lsr = −10 km s −1 . To avoid contribution from these background clouds, the emission between −6.5 km s −1 to −4 km s −1 and −11 km s −1 to −9 km s −1 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 2 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 12 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.
Notes. The last three columns display the total outflow momentum rate, assuming symmetric outflow, corrected for the CO opacity at an outflow inclination of 57.6 • , 70 • and 80 • . The last two columns are only given for the bipolar outflow originating in the H1 fragment as it is expected that it might have such a large inclination angle with the line of sight.
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. On 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 12 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 12 CO(2-1) data.

The outflow momentum rate
Integrating the 12 CO(2-1) wing emission over the entire extent of the red-and blue-shifted outflows, one can estimate their CO outflow momentum rate assuming the driving source is located at a velocity of −21 km s −1 . 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 redshifted 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 . 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 −3 M km s −1 yr −1 , 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 or 4σ C 18 O (instead of 3.5σ C 18 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 CO 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 −3 and 3.1 × 10 −2 M km s −1 yr −1 , 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.

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 12 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 12 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 12 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 • .

H II emission in the cavities
In the previous section we showed that the presence of midinfrared 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  provides radio continuum data at 4.85 GHz (6.19 cm) at an angular resolution of 5 (Condon et al. 1993). In this 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 −1 ).
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. Only 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 H II 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 −1 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 −1 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 2 cm −3 and n e < 4.6 × 10 2 cm −3 , 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 6 and 6.9 × 10 6 K cm −3 in N1 and N2, respectively. Compared to the thermal and ram pressure of the molecular gas, estimated to be ∼2.6 × 10 5 K cm −3 (assuming T K = 20 K) and 1.2 × 10 7 K cm −3 (assuming v = 2 km s 1 ) 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 2 cm −3 and resulting thermal H II pressure <3.7 × 10 6 K cm −3 .
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;Rodríguez-Kamenetzky et al. 2017). Dedicated radio continuum observations are required to resolve this question.

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.

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 selfconsistent 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. 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).

Dust continuum images
When the photon packages leave the volume grid, they produce an image at different wavelengths in units of Jy pixel −1 . For each of the simulations 10 7 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 −1 which can be easily converted to MJy sr −1 units and compared with the observations.

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).

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) where n c is the central density. The power law index η was taken to be 2 (Schmiedeke et al. 2016), and | r | is given by 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 with: In these equations, R is the radial distance in the equatorial plane, n 0 is the characteristic density, r 0 the characteristic radius, β the disk flaring power, q the surface density radial exponent, and h 0 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 5 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 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.

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 A22, page 10 of 22 L. Bonne et al.: A potential new phase of massive star formation. A luminous outflow cavity centred on an infrared quiet core

Model
Plummer core Passive flared disc (PFD) Cavities Model

Notes.
A combination of all of these input parameters were run. In Model 3, r 0,x of the Plummer core and r 0 of the PFD are fixed to be the same.
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.

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 A22, page 11 of 22 A&A 665, A22 (2022)  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 4 AU, a characteristic density of 10 7 cm −3 and a cavity density of 10 4 cm −3 . This is to be directly compared to Fig. 3 (left). Note the order of magnitude difference in flux density between the two images.
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.

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 χ 2 calculation: where S core λ are the core fluxes estimated from the synthetic images at the wavelength λ (i.e. 250, 350 and 500 µm, respectively). O 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 χ 2 core values as a function of the central source luminosity. Only models achieving χ 2 core < 3 are displayed. At first sight, there appears to be no clearly preferred luminosity, both 500 L and 4000 L runs achieve low χ 2 core values. However, at 4000 L these low χ 2 core 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 χ 2 core ≤ 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.

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.

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 [OI] 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(Lehmann et al. , 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 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(Lehmann et al. , 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 CO entrainment efficiency of 50% (i.e. f ent ; Duarte-Cabral et al. 2013) and the observed CO momentum rate at an 80 • inclination angle (F CO ). This gives a jet velocity of 3.8 × 10 3 km s −1 . This is significantly above the typical velocity range of 100-1000 km s −1 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 18 O(2-1) and 13 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 13 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 18 O(2-1) and 13 CO(2-1) cover up to 8-10 km s −1 . 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 13 CO(2-1) and C 18 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-ofsight velocity gradients. High-angular resolution observations of well-known shock tracers will be necessary to resolve this issue.

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 protostellar 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 12 CO(2-1), 13 CO(2-1), and C 18 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 13 CO(2-1) spectrum, which gives a derived infall velocity of 0.5 km s −1 (using 12 CO(2-1) gives a similar value of 0.6 km s −1 ). 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 −1 . 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 18 O(2-1) hub centroid velocity obtained in C 18 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. 2007Peretto et al. , 2013Smith 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 −1 pc −1 , while the average value over the hub is 1.9 +0.3 −0.2 km s −1 pc −1 . 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. 2013Henshaw et al. , 2014Peretto 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 freefall time. The inflow timescale is given by: where R is the hub radius. Taking R = 1 pc and v inf = 1 km s −1 , we obtain t inf = 1 Myr. On the other hand, the free fall time scale is given by: 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 −20 g cm −3 which corresponds to a molecule number density of n H 2 = 1.2 × 10 4 cm −3 . 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: 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 ρ is reached at a radius that is about half the external hub radius of 1 pc. With these values in hand we obtainṀ hub inf = 2.2 × 10 −3 M yr −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 −1 with apparent peaks as high as 250 • kyr −1 . 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 3 and 2.7 × 10 4 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.

Summary and future prospects
We have presented the first study of the G345.88-1.10 hub filament system (d = 2.26 +0.30 −0.21 kpc, M = 2.9 × 10 3 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 12 CO(2-1), 13 CO(2-1), C 18 O(2-1) and H30α towards the hub. The main results are: -The H1 fragment has an estimated mass of 2.1 +1.2 −1.1 × 10 2 M within a radius R eff = 0.14 pc, and a luminosity upper limit of 10 2 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 3 L . It is believed that the main axis of the bipolar nebulosity is lying close to the plane of the sky.
-12 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 −2 M km s −1 yr −1 . The nebulosities thus trace the outflow cavity of a powerful protostellar outflow emanating from the H1 fragment with an opening angle of ∼90 ± 15 • . A22, page 15 of 22 A&A 665, A22 (2022) -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 −3 M yr −1 ).
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 ; iii. map the spatial distribution of shocked gas via the detection of SiO 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, SOFIA observations of [C II], [O I], H 2 , 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. 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 χ 2 core (see Fig. G.1). First, we find that only a small number of runs with nebulosity have a χ 2 core < 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 L , the central density n c of the Plummer profile, or the characteristic density n 0 of the passive flared disc, needs to be larger than n H 2 = 5×10 6 cm −3 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 ≥ 10 7 cm −3 . 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.

Appendix H: Calculating the cooling timescale for the H1 fragment and the cavities
The cooling time scale is given by 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 = 3 2 k b T gas µ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 A22, page 21 of 22 A&A 665, A22 (2022) T gas = 20 K T gas = 40 K T gas = 60 K T gas = 80 K T gas = 100 K T dust = 20 K 4.4 yr 8.8 yr 13 yr 18 yr 22 yr T dust = 30 K 3.9×10 −1 yr 7.7×10 −1 yr 1.2 yr 1.5 yr 1.9 yr T dust = 40 K 6.9×10 −2 yr 1.4×10 −1 yr 2.1×10 −1 yr 2.7×10 −1 yr 3.4×10 −1 yr proton. The last term that has to be constrained is the cooling rate per unit volume (λ cool ). If the gas is dense enough, typically n H 2 > 10 5 cm −3 (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×10 5 cm −3 , 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 where κ ν is the dust opacity, D the dust to gas mass fraction and B ν the Planck function. Using Ossenkopf & Henning (1994) this results in λ cool = 4.68 × 10 −31 nT 6 dust ( in erg s −1 cm −3 ) (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. 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×10 4 cm −3 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 cm −3 , λ cool is given by λ cool = 1.5 × 10 −26 T 2.7 ( in erg s −1 cm −3 ) (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 λ cool between 4.9×10 −23 -3.8×10 −21 erg cm −3 s −1 which is 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 3 and 2.7×10 4 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.