Dense gas formation in the Musca filament due to the dissipation of a supersonic converging flow

Observations with the Herschel Space Telescope have established that most of the star forming gas is organised in interstellar filaments, a finding that is supported by numerical simulations of the supersonic interstellar medium (ISM) where dense filamentary structures are ubiquitous. We aim to understand the formation of these dense structures by performing observations covering the $^{12}$CO(4-3), $^{12}$CO(3-2), and various CO(2-1) isotopologue lines of the Musca filament, using the APEX telescope. The observed CO intensities and line ratios cannot be explained by PDR (photodissociation region) emission because of the low ambient far-UV field that is strongly constrained by the non-detections of the [C II] line at 158 $\mu$m and the [O I] line at 63 $\mu$m, observed with the upGREAT receiver on SOFIA, as well as a weak [C I] 609 $\mu$m line detected with APEX. We propose that the observations are consistent with a scenario in which shock excitation gives rise to warm and dense gas close to the highest column density regions in the Musca filament. Using shock models, we find that the CO observations can be consistent with excitation by J-type low-velocity shocks. A qualitative comparison of the observed CO spectra with synthetic observations of dynamic filament formation simulations shows a good agreement with the signature of a filament accretion shock that forms a cold and dense filament from a converging flow. The Musca filament is thus found to be dense molecular post-shock gas. Filament accretion shocks that dissipate the supersonic kinetic energy of converging flows in the ISM may thus play a prominent role in the evolution of cold and dense filamentary structures.


Introduction
Observations with the Herschel Space Telescope have revealed that filamentary structures are ubiquitous in the supersonic interstellar medium (ISM) (e.g. André et al. 2010;Molinari et al. 2010;Henning et al. 2010;Arzoumanian et al. 2011;Schneider et al. 2012). However, there is an ongoing discussion regarding the nature and diversity of the filaments, namely, whether they are sheets viewed edge-on or, rather, dense gas cylinders. It is also considered whether we observe a full range of filament classes: from cross-sections of sheets to dense, star-forming cylindrical structures. In any case, understanding the nature, formation and evolution of filaments is essential as they are the sites of star formation. This was demonstrated by a number of recent studies which showed that pre-and protostellar cores are mostly located in filaments (e.g. Polychroni et al. 2013;André et al. 2014;Schisano et al. 2014;Könyves et al. 2015;Marsh et al. 2016;Rayner et al. 2017). In numerical simulations of the ISM, filaments are omnipresent and can form in various ways: through shocks in (magnetic) supersonic turbulent colliding flows (e.g. Padoan et al. 2001;Jappsen et al. 2005;Smith et al. 2016;Federrath 2016;Clarke et al. 2017;Inoue et al. 2018), during the pc (Franco 1991). We adopt the nomenclature from Cox et al. (2016) for the different features of the Musca filament and we indicate them in Fig. 2 (see also Fig. 2 in Cox et al. (2016)). First, there is the high column density filament crest with N>3 10 21 cm −2 (similar value to the one used by Cox et al. (2016) which is N>2.7 10 21 cm −2 ). In addition, C 18 O(2-1) emission basically disappears below this column density value. Attached to the filament crest are intermediate column density (N∼2 10 21 cm −2 ) hair-like structures called strands, with a size of ∼0.2-0.4 pc. Even further outwards at larger distances, narrow straight structures of up to a few parsec lengths, called striations (Goldsmith et al. 2008;Palmeirim et al. 2013;Alves de Oliveira et al. 2014;Heyer et al. 2016;Cox et al. 2016;Malinen et al. 2016) are seen. They are mostly orthogonal to the crest and located in the ambient cloud. The ambient cloud then refers to the environmental gas, traced by the extinction map in Fig. 1, embedding the denser Musca filament crest and strands. The filament crest is velocity-coherent ) and has only one young stellar object (YSO) located at the northern end (e.g. Vilas-Boas et al. 1994, Fig. 1). It was recently proposed that Musca is a sheet seen edge-on (Tritsis & Tassis 2018). However, one of the main results of a companion paper (Bonne et al., in prep., hereafter Paper I), is that the Musca filament crest is more consistent with a cylindrical geometry at a density of n H 2 ∼ 10 4 cm −3 . Because other authors (Cox et al. 2016;Kainulainen et al. 2016;Hacar et al. 2016) also support this cylindrical geometry, we adopt this view as a working hypothesis.
In Paper I, the large scale kinematics and physical conditions in the Musca filament and cloud are studied with the main CO(2-1) isotopologues, as well as the relation between Musca and the Chamaeleon-Musca complex. In this paper, we present 12 CO(4-3) and 12 CO(3-2) lines and argue for the presence of a blueshifted excess velocity component best visible towards the filament crest and strands. A detailed analysis of this blueshifted component is the objective of this paper. Here, we present a non-LTE CO excitation analysis of the blueshifted component, as well as complementary observations of the far-infrared fine structure transitions of carbon ([C I]), oxygen ([O I]) and ionised carbon ([C II]), which reveals a warm gas component around the filament. We show that the excitation of this warm gas can be explained by a filament accretion shock, that is, a low-velocity shock due to mass accretion on a filament from a converging flow in an interstellar cloud.

APEX
We used the FLASH460 receiver on the APEX telescope (Güsten et al. 2006;Klein et al. 2014), and obtained a single pointing in [C I] 3 P 1 -3 P 0 at 609 µm towards the Musca filament at α 2000 = 12 h 28 m 55 s and δ 2000 = -71 o 16 55 , using the OFF position α 2000 = 12 h 41 m 38 s and δ 2000 = -71 o 11 00 ). The beamsize is 13 and the spectral resolution is ∼ 0.05 km s −1 . The data reduction was done with the CLASS software 1 . A main beam efficiency 2 η mb = 0.49 was applied for [C I], a baseline of order one was removed, and a correction for a 490 kHz shift in the FLASH460 instrument was performed (see Paper I, F. Wyrowski, priv. comm.). Fitting the [C I] single pointing (Fig. 3) with a Gaussian profile provides a peak temperature brightness T mb = 3.5 K and FWHM = 0.8 km s −1 (or 3.4·10 −7 erg s −1 cm −2 1 http://www.iram.fr/IRAMFR/GILDAS 2 http://www.apex-telescope.org/telescope/efficiency/ sr −1 ). The CO observations with APEX are presented in more detail in Paper I, but here we shortly summarise the most important facts. All data was obtained in 2017 and 2018, using FLASH345 and FLASH460 for the 12 CO(3-2) and 12 CO(4-3) mapping of the northern and southern regions (Fig. 1). The FLASH345 (FLASH460) observations have a spectral resolution of 0.033 (0.05) km s −1 and an angular resolution of 18 (14 ). Main beam efficiencies of η mb =0.65 (0.49) were applied to the antenna temperatures. We here use data sampled to 0.1 km s −1 and smoothed to a resolution of 28 . This allows a better comparison to the CO(2-1) APEX data which were taken with the PI230 receiver. The spectral resolution for CO(2-1) is 0.08 km s −1 at an angular resolution of 28 . Here, we applied a main beam efficiency of η mb =0.68. The center (0,0) position of the maps is α 2000 = 12 h 28 m 58 s , δ 2000 = -71 • 16 55 for the northern map and α 2000 = 12 h 24 m 46 s , δ 2000 = -71 • 47 20 for the southern map.

SOFIA
In June 2018, the upGREAT instrument (Risacher et al. 2018) on the Stratopsheric Observatory for Far-Infrared Astronomy (SOFIA) was used for a 70 minute single pointing of the 7 pixel array covering the atomic oxygen [O I] 3 P 1 -3 P 2 fine structure line at 63 µm and the ionised carbon [C II] 2 P 3/2 -2 P 1/2 fine structure line at 158 µm towards the southern APEX map at α 2000 = 12 h 24 m 41 s .6, δ 2000 = -71 • 46 41 . 0. These observations were performed in the single beam switching mode with a chop amplitude of 150 towards a location with weak or no emission in all Herschel far-infrared bands (70-500 µm). This was done to obtain a good baseline for such sensitive observations. Mars was used as a calibrator to determine the main beam efficiencies of the individual pixels. The fitted water vapor column was typically around 10 µm. All the intensities reported here are on the main beam temperature scale. The [O I] and [C II] observations have a main beam size of 6.3 and 14.1 , respectively, and both data sets were smoothed to a spectral resolution of 0.3 km s −1 . There is no detection of the [C II] and [O I] lines (Sec. A) in the individual pixels ( Fig. A.1) or in the array-averaged spectra (Fig. 3). The 3σ upper limits for the two lines are 0.11 K and 0.19 K, equivalent to 7.5·10 −7 erg s −1 cm −2 sr −1 and 2.1·10 −5 erg s −1 cm −2 sr −1 when assuming a FWHM of 1 km s −1 , for [C II] and [O I], respectively.

Results and analysis
3.1. The ambient FUV field from a census of nearby stars We estimated the FUV field upper limit in Musca using a census of nearby ionising stars (e.g. Schneider et al. 2016). For that, we employed the SIMBAD 3 database (Wenger et al. 2000) and imposed a search for B3 and earlier type stars within a radius of 40 • (100 pc), centered on the Musca filament, and a maximal distance of 250 pc. This provided a census of 59 ionising stars including two dominant O9-type stars. Knowing the spectral type of these ionising stars, we estimated the contribution of the FUV emission using model atmospheres 4 of Kurucz (1979) with solar abundances (Grevesse & Sauval 1998) and log(g) ∼4.0 (g is the surface gravity). The temperature of the different spectral types was obtained from Fitzpatrick & Massa (2005) and Pecaut & Mamajek (2013). The distance information and the location  (Cox et al. 2016) inserted into the extinction map of the Musca cloud (Schneider et al. 2011), which is scaled to the Herschel column density map. The filament crest is defined as N > 3·10 21 cm −2 (purple and white). The region outlined by a red rectangle is the zoom displayed in Fig. 2. The ambient cloud is displayed in green and blue. The northern and southern areas mapped with FLASH+ on the APEX telescope are indicated in black. The triangle indicates the location of the only YSO in the Musca filament. on the sky of both Musca and the stars then allows to calculate their relative distances. We assumed a r −2 decrease of the flux and projected all stars in the plane of the sky, and obtained in this way an ambient FUV field of 3.4 G 0 (Habing 1968) for the Musca cloud as an upper limit. The real field is lower because we did not take into account any extinction, that is, attenuation by diffuse gas and blocking by molecular clumps.

Limits on FUV heating from CI & CII
As a second method to constrain the ambient FUV field, we used the SOFIA [C II] upper limit and the APEX [C I] brightness observed in Musca. [C II] is an excellent tracer of the FUV field at low densities as it is a direct result of the ionising radiation (e.g. Tielens & Hollenbach 1985) and its non-detection points towards a low value for the ambient FUV field. [C I] is not a direct product of FUV ionisation, but PDR models demonstrate that its brightness is also sensitive to the FUV field (e.g. Hollenbach & Tielens 1997;Röllig et al. 2007). In this section, we compare the observed line brightness with predictions from PDR models using the 2006 version of the PDR toolbox 5 (Kaufman et al. 2006;Pound & Wolfire 2008). These plots are expressed in molecular hydrogen density n H 2 , where we assume that all hydrogen atoms are locked in molecular hydrogen.
The [C II] upper limit restricts the FUV field to values < 1 G 0 for densities n H 2 ≤ 10 4 cm −3 . This is consistent with the intensity of the [C I] line which is so weak that it is below the minimum value in the PDR toolbox. Even with a beam filling value of 30 % for [C I], which would be low around the Musca filament, the FUV field remains below 1 G 0 for typical densities of the Musca ambient cloud. Using the Meudon PDR code 6 (Le Petit et al. 2006;Le Bourlot et al. 2012;Bron et al. 2014) for the same lines, confirms that the upper limit on the [C II] intensity can only be the result of a FUV field strength < 1 G 0 . With these observational upper limits for the ambient FUV field, the A&A proofs: manuscript no. template  Kaufman et al. (2006) PDR models restrict the maximal surface temperature of the cloud to 25 K. This temperature upper limit for gas embedded in such a weak FUV field is also found in other theoretical models where heating by the FUV field is taken into account (e.g. Godard et al. 2019), and consistent with the maximal CO temperature in numerical simulations of molecular clouds evolving in a FUV field ≤ 1.7 G 0 Clark et al. 2019). The spectra of the CO isotopologues unveil three velocity components that show small velocity variations in the north and south. First, there is a single component in C 18 O(2-1) at 3.5 km s −1 (north) and 3.1 km s −1 (south) that traces the Musca filament crest . In 13 CO(2-1) one observes a blueshifted shoulder to the velocity component of the filament crest, which is a velocity component related to the strands, namely, the dense interface region between the filament crest and the ambient cloud. It has a typical velocity of 3.1 km s −1 in the north and 2.7 km s −1 in the south and is not detected in C 18 O with current data. In 12 CO a third velocity component is observed, without any clearly detected counterpart in 13 CO or C 18 O, which is further blueshifted so that we call it the blueshifted component . The brightness peak of the blueshifted component in the northern map occurs at 2.8 km s −1 and in the southern map at 2.5 km s −1 . This blueshifted component is the focus of this paper.

Three CO velocity components
The individual 12 CO(4-3) spectra in the northern map are presented in Fig. 4, overlaid on the Herschel column density and temperature map in the top panels and together with individual spectra in other CO isotopologues at selected positions in the bottom panel. Inspecting the 12 CO(3-2) and 12 CO(4-3) spectra, we find that these lines are clearly detected towards the lower column density regions and that this blueshifted component is present over the entire map, see Fig. 4.
Because our assumption is that the Musca filament is cylindrical, the observed CO emission in direction of the crest and strands also contains gas from the ambient cloud along the lineof-sight. In particular, the 12 CO(2-1) blueshifted line can have a more significant component arising from this gas phase. We can thus not fully exclude effects of high optical depth and selfabsorption, so we run tests using a two-layer gas model for the CO(2-1) isotopologue lines in order to calculate a possible impact of foreground absorption and present the results in Appendix B. For that, we assumed only 2 CO line components, i.e. the crest component and a shoulder+blueshifted one that it self-absorbed and only apparently shows separate components. However, the result of this modelling is that it is not possible to reproduce the observed CO line profiles with only two components, so that we are confident in our approach to use 3 separate lines. The best fitting model for the blueshifted component of all CO lines in the northern and southern map is the one with the parameters given in Table 1. We note from the Table that the fitted 12 CO(2-1) linewidth of the blueshifted component in the south is higher than the linewidth obtained for 12 CO(3-2) and 12 CO(4-3), which have a similar width. It is thus possible that there, the low-J 12 CO(2-1) line traces more material along the line-of-sight from the ambient cloud. Inspecting the spectra, see Fig. 3, confirms that 12 CO(3-2) and 12 CO(4-3) only show up at higher velocities and have a smaller linewidth. This difference, and the irregular 12 CO(2-1) shape in the southern map, suggests that 12 CO(2-1) traces the ambient cloud down to lower velocities than the detected 12 CO(4-3) and (3-2) emission. We also determined from the averaged 12 CO(2-1) and 13 CO(2-1) spectra that the 12 CO(2-1)/ 13 CO(2-1) line ratio in the blueshifted component in the northern and southern map are between 15 and 60, which are expected values for optically thin emission. Note that fractionation (see Appendix B) can cause observed values for the 12 CO/ 13 CO ratio significantly below 60. This is in agreement with early indications for Musca of significant CO isotopologue abundace variations ) which we further confirm in Paper I.

Density for the blueshifted component
From the sections before, we learned that the blueshifted component is an individual feature towards the crest and strand regions, and that the 12 CO(2-1) line has contributions from the ambient cloud along the line-of-sight. Here, we estimate the typical density and upper limit for this ambient cloud close to the filament. For that we use the density of the fitted Plummer profile at the outer radius of the filament. With the values from Cox et al. (2016), correcting for a distance of 140 pc for the Musca filament, this gives n H 2 = 4.0·10 2 cm −3 at r = 0.2 pc. An approach to estimate the density upper limit is to combine the maximal column density associated with the ambient cloud, which is N∼10 21 cm −2 (Cox et al. 2016), with the minimal possible size of the ambient cloud. As it is unlikely that the ambient cloud has a smaller size along the line of sight than the filament, this gives a minimal size of 0.4 pc. Using a minimal size of 0.4 pc along the line of sight gives an upper limit of n H 2 = 8·10 2 cm −3 for the ambient cloud.
The [C I] emission at the velocity of the blueshifted component in Fig. 3 is weak, ∼1-1.5 K, but it is difficult to constrain this because of the noise. This indicates that [C I] can trace the ambient cloud down to at least n H 2 ∼ 4·10 2 cm −3 in a weak FUVfield. More extensive [C I] observations will be required to better understand the physical conditions traced by [C I] in the ISM. This is particularly important for comparison with simulations since the physical conditions traced by [C I] remain uncertain in simulations (e.g. Glover et al. 2015;Franeck et al. 2018;Clark et al. 2019).

12 CO(4-3) excess emission in the blueshifted velocity component
The 12 CO(3-2) and 12 CO(4-3) transitions have relatively high excitation temperatures (33 K and 55 K respectively; Müller et al. 2005). This allows us to investigate the presence of warm gas in the Musca cloud. We focus on the blueshifted component for which 12 CO(4-3) is observed over the full map. This could be an indication of quite uniform heating by e.g. a FUV field, though there are indications of spatial variations in the data. However, because of the noise in the spatially resolved map, we will not focus on these possible spatial variations in this paper. To further investigate the heating that leads to the 12 CO(4-3) emission, we study the 12 CO(4-3)/ 12 CO(2-1) and 12 CO(4-3)/ 12 CO(3-2) brightness temperature ratios. In the previous section, we noted that it is important to take care comparing the integrated brightness of 12 CO(2-1) with 12 CO(3-2) and 12 CO(4-3) in the blueshifted component because of their difference in linewidth. Furthermore, the noise in the 12 CO(3-2) and 12 CO(4-3) spectra makes it impossible to confidently fit a spectrum for every pixel in the 12 CO(4-3)/ 12 CO(2-1) 12 CO(4-3)/ 12 CO(3-2) mean median mean median north (2.8 km s −1 ) 0.38 ± 0.06 0.38 0.62 ± 0.09 0.61 south (2.5 km s −1 ) 0.39 ± 0.13 0.38 0.55 ± 0.19 0.53 Table 2. Mean and median main beam brightness ratios for 12 CO(4-3)/ 12 CO(2-1) and 12 CO(4-3)/ 12 CO(3-2) in the blueshifted component for the northern and southern map.
map, even when spatially smoothing the data. Because of this, we will approach this section by focusing on the peak brightness temperatures.
The peak brightness of the blueshifted component occurs near 2.8 km s −1 in the northern map and near 2.5 km s −1 in the southern map. At these velocities, we determined the brightness temperature ratios for every pixel in the maps. In Fig. 5, the distribution of the 12 CO(4-3)/ 12 CO(2-1) and 12 CO(4-3)/ 12 CO(3-2) peak brightness ratios for the northern map are displayed. We determined the average and median values for the various ratios in both maps which are summarized in Tab. 2. Overall, there is no large difference between average and median values. For the 12 CO(4-3)/ 12 CO(2-1) and 12 CO(4-3)/ 12 CO(3-2) ratio in the northern map, we obtain a median value of 0.38 and 0.61 at 2.8 km s −1 , respectively. In the southern map, the noise of both the 12 CO(4-3) and 12 CO(3-2) transition is higher as it was observed under poorer weather conditions. The median 12 CO(4-3)/ 12 CO(2-1) ratio is 0.38 at 2.5 km s −1 and the respective value for the 12 CO(4-3)/ 12 CO(3-2) ratio is 0.53.

Modeling with the PDR Toolbox
We compared the observed 12 CO(4-3)/ 12 CO(2-1) and 12 CO(4-3)/ 12 CO(3-2) brightness temperature ratios in the blueshifted component with predictions using the PDR toolbox in the allowed density range of the ambient cloud in Musca, see Sec. 3.3.2. Figure 6 shows that for a FUV field strength <1 G 0 , the predicted 12 CO(4-3)/ 12 CO(2-1) brightness temperature ratio is smaller than 0.05, which is more than a factor 5 lower than observed in Musca. The same is found for the 12 CO(4-3)/ 12 CO(3-2) brightness temperature ratio, with predicted ratios of ∼ 0.1, Fig. 6. Observed 12 CO(4-3)/ 12 CO(2-1) and 12 CO(4-3)/ 12 CO(3-2) ratios compared with predictions by the PDR Toolbox. The shaded areas, constructed by logarithmic interpolation between the points in the PDR Toolbox density grid, show the evolution of the 12 CO(4-3)/ 12 CO(2-1) and 12 CO(4-3)/ 12 CO(3-2) brightness temperature ratio as a function of the FUV field in the allowed density range of the Musca cloud. At a typical density n H 2 = 4·10 2 cm −3 of the ambient cloud (dashed line) for the allowed FUV field strength (< 1 G 0 ), the predicted line ratios are only a fraction of the observed ones in the blueshifted component.
which is again more than a factor 5 lower than observed. It is thus impossible to obtain the observed ratios for the densities and FUV field strength in the ambient Musca cloud that are allowed by [C II], [C I], [O I], and the calculated upper limit from the nearby census of OB stars. This strongly suggests that the 12 CO(4-3) emission in the low-column density blueshifted component cannot be explained as a result of FUV heating.

Modeling with RADEX
We further investigate whether the 12 CO(4-3)/ 12 CO(2-1) and 12 CO(4-3)/ 12 CO(3-2) brightness temperature ratios in the blueshifted component can be reproduced with the non-LTE RADEX code (van der Tak et al. 2007). We use a non-LTE approach because the high critical density (n H 2 > 5·10 4 cm −3 ) of 12 CO(3-2) and 12 CO(4-3) implies that these lines are subthermally excited. For RADEX analysis we use a FWHM of 0.4 km s −1 , temperatures between 15 and 25 K, and representative densities for the ambient cloud: n H 2 = 3·10 2 , 5·10 2 , 7.5·10 2 and 1.5·10 3 cm −3 . For each RADEX model, we also need an upper limit on the 12 CO column density of the blueshifted component. This can be done by calculating with RADEX the 13 CO column density from the 13 CO(2-1) brightness (∼0.6 K, see Fig. 3) for every density and temperature. In a weak FUV field the [ 12 CO]/[ 13 CO] abundance ratio is ≤ 60 (Visser et al. 2009;Röllig & Ossenkopf 2013), which puts an upper limit on the 12 CO column density and opacity. For further analysis, we use RADEX models with a predicted 12 CO(2-1) brightness up to 50% brighter than observed towards Musca to take into account some uncertainties such as possible opacity broadening of 12 CO(2-1) compared to 12 CO(4-3), a non unity beam filling or different line calibration (uncertainty in η mb ). Studying the 12 CO(4-3)/ 12 CO(2-1) and 12 CO(4-3)/ 12 CO(3-2) brightness temperature ratios within these limits, we find that the ratio increases with increasing density, column density Fig. 7. 12 CO(4-3)/ 12 CO(2-1) brightness temperature ratios predicted with RADEX for the plausible range of physical conditions in the ambient cloud of the Musca filament. The lines connect the RADEX models with identical density and temperature for varying 12 CO column densities. The predicted ratios with RADEX show a huge discrepancy with the observed 12 CO(4-3)/ 12 CO(2-1) ratio towards the Musca cloud (indicated by the dashed line). The standard deviation of the observed values from the average is indicated by the grey area. The same is observed for the 12 CO(4-3)/ 12 CO(3-2) brightness temperature ratio. or temperature, see Fig. 7. However, no model is capable to reproduce more than 40% of the observed 12 CO(4-3)/ 12 CO(2-1) and 12 CO(4-3)/ 12 CO(3-2) ratios towards the Musca ambient cloud, see Fig. 7 for the 12 CO(4-3)/ 12 CO(2-1) ratio. Taking higher densities for the ambient cloud than put forward in sect. 3.3.2, does not offer a solution for bringing calculations in agreement with observations. Though this can increase the predicted ratio for a fixed temperature and 12 CO column density, this also strongly increases the line brightness temperature. Consequently, one either has to reduce the temperature or 12 CO column density to keep line brightnesses that are not too far off from the observed values. This effort again forces brightness temperature ratios to values found in Fig. 7.

The CO line ratio with RADEX
In order to increase the 12 CO(4-3)/ 12 CO(2-1) and 12 CO(4-3)/ 12 CO(3-2) brightness temperature ratio while keeping the brightness of both lines low enough, one needs to consider a small layer embedded in the diffuse gas that can increase the 12 CO(4-3) brightness without significantly increasing the 12 CO(2-1) and 12 CO(3-2) emission. Based on the 12 CO(4-3) excitation conditions, the layer should contain warm gas (> 50 K). In fact, this 'layer' can be clumpy, that is, pockets of warm gas that are embedded in more tenuous interclump gas. However, modelling such a scenario is not possible with RADEX and out of the scope of this paper. Our objective here is to show the existence of a warm gas component. We thus investigated with RADEX the impact of a warm gas layer by running models with A&A proofs: manuscript no. template Fig. 8. Predicted brightness temperature ratio of 12 CO(4-3)/ 12 CO(2-1) by RADEX as a function of temperature for different densities. This demonstrates the need for a warm and dense CO layer to at least reach the observed brightness ratios in the blueshifted component. For each density we plot the ratio for a column density such that the predicted brightness of 12 CO(4-3) by RADEX is roughly similar to the observed brightness towards Musca (T mb =1.5 -3 K) at the temperatures that obtain sufficiently high brightness ratios. The lowest three densities (n H 2 = 10 3 , 3·10 3 and 5·10 3 cm −3 ) use N12 CO = 1.1·10 15 cm −2 , and n H 2 = 7·10 3 cm −3 uses N12 CO = 9·10 14 cm −2 . The black horizontal line indicates the average observed ratio towards Musca, and the grey area indicates the standard deviation. Fig. 9. Predicted brightness temperature ratio of 12 CO(4-3)/ 12 CO(3-2) by RADEX for the same models that are presented in Fig. 8. It confirms the need for high temperatures (> 50 K) and densities (5-7·10 3 cm −3 ) to reach the observed excitation of 12 CO(4-3). The black horizontal line indicates the average observed ratio towards Musca, and the grey area indicates the standard deviation. a FWHM of 0.4 km s −1 as an overall average, and temperatures between 10 and 150 K at different densities. Note that a temperature of 150 K would result in a thermal FWHM of 0.5 km s −1 . This is higher than the observed typical FWHM of 0.4 km s −1 , but we anticipate that there is also uncertainty on the FWHM so that we work with a kinetic temperature upper limit of 150 K. We additionally restrain ourself to models that manage to reproduce the observed 12 CO(4-3) brightness temperature in Musca 7 . Figures 8 and 9 show the results for different densities. Both the 12 CO(4-3)/ 12 CO(2-1) and 12 CO(4-3)/ 12 CO(3-2) ratios indicate that a temperature > 50 K as well as relatively high densities are required to reach the observed ratios. The 12 CO(4-3)/ 12 CO(3-2) ratio, in particular, points to high densities and temperature, see Fig. 9, in order to reach the observed ratios.
Though the RADEX analysis demonstrates the need for warm and dense gas to obtain the observed 12 CO(4-3) brightness in the blueshifted component, we can not narrow down more precisely the temperature and density range because there are no higher-J CO observations at the moment. This warm and dense gas exactly fits with the predictions for gas heated by lowvelocity shocks (Pon et al. 2012;Lesaffre et al. 2013), and socalled slow-type (with regard to their phase velocity) magnetised shocks with v s = 1-3 km s −1 can easily reach temperatures as high as 100 K . These slowtype shock models predict a physical size of the warm gas layer around 10 15 cm, which fits with the estimated physical sizes for the RADEX models. The estimated size, that is the layer thickness, of the models that manage to reach the observed brightness temperature ratios in Figs. 8 and 9 are 1.3·10 15 cm and 2.2·10 15 cm, respectively at n H 2 = 7·10 3 cm −3 (with N12 CO = 9·10 14 cm −2 ) and 5·10 3 cm −3 (with N12 CO = 1.1·10 15 cm −2 ). To calculate the physical size of the RADEX models an abundance of H 2 / 12 CO = 10 4 was used, which is a typical value for weakly irradiated molecular gas. Lastly, we note that this warm gas is observed in both maps and thus likely universally present around the Musca filament.

CI column density in the blueshifted component
We noted in Fig. 3 the presence of [C I] emission in the blueshifted component with a brightness temperature of the order of 1-1.5 K. From this weak [C I] emission we here estimate the [C I] column density. We use a FWHM = 0.4 km s −1 , a temperature of 20 K, and typical densities n H 2 = 5-7.5·10 2 cm −3 for the ambient Musca cloud with RADEX. This points to N CI ∼ 10 16 cm −2 , see Tab. 3. Comparing this with the same models for 12 CO emission from the ambient cloud in section 3.4, we find that at least 20% and possibly up to 50% of carbon is still found in its atomic form.
On the other hand, it is observed from Tab. 3 that the warm and dense gas layer, necessary to explain the bright 12 CO(4-3) emission, provides a negligible contribution to the [C I] emission. We thus emphasise again that the blueshifted component has two contributions: -The ambient cloud which gives rise to [C I] and low-J CO emission.
-A warm gas layer (or pockets of warm gas) with little [C I] emission that is responsible for the bright CO(4-3) line.

Shock models
For an in-depth comparison with shock models, it would be preferable to have additional observations of mid-to high-J (J up > 4) CO lines (e.g. Pon et al. 2012 2016) since 12 CO(2-1), and possibly some 12 CO(3-2), emission in the blueshifted component also comes from non-shocked gas (as it is found in synthetic observations of simulations; see , it is not well investigated to which extent mid-J CO lines, in particular the 12 CO(4-3) transition, contribute as a cooling line for shocks. We thus compared our observations to results of the Paris-Durham shock code (Flower & Pineau des Forêts 2003;Lesaffre et al. 2013;Godard et al. 2019) that are available for both non-irradiated C-and J-type shocks. Looking into computed model grids, we find that non-irradiated C-type shocks with a pre-shock density n H 2,0 ∼ 5·10 2 cm −3 do not manage to reproduce the observations. Both the integrated intensity of 12 CO(4-3) and the predicted ratios are below the observed values in Musca. However, J-type shocks, which are good first order models for slow-type shocks , with a pre-shock density n H 2,0 ∼ 5·10 2 cm −3 fit with the 12 CO observations, see Fig. 10. The predictions by these models for the [C II], [O I], and [C I] brightnesses are also in agreement with the observations of Musca (see Appendix A). J-type shocks as a proxy for slow-type low-velocity shocks are justified since the dynamics of these slow-type shocks are driven by the gas pressure. Note that we used published shock models from the Paris-Durham code for this first comparison which have shock velocities of 4 km s −1 and higher, while in Musca the shock velocity can be lower. A comparison with shock models at lower velocities is work in progress, and will be addressed in a forthcoming paper. The nondimensional parameter b = (B/1µG)/ √ n cm −3 = 0.1 in the shock models, with B the component of the magnetic field that is perpendicular to the direction of shock propogation, is used to model J-type shocks, and b = 1 is used to model C-type shocks (Lesaffre et al. 2013).
To perform a first estimate of the shock velocity, we have to take into account that there is an orientation angle for the shock propagation with respect to the plane of the sky (POS). We will assume shock propagation along the magnetic field. In this scenario, one can estimate the true shock velocity when the angle between the magnetic field and the POS is known, using where γ is the angle between the magnetic field and the POS, and v los is the observed shock velocity along the line of sight for which we take 0.4 km s −1 . In Planck Collaboration et al. (2016), Fig. 10. 12 CO(4-3)/ 12 CO(2-1) integrated brightness ratio (black) for Jtype shock models from the Paris-Durham code. The ratio at a shock velocity of 4 km s −1 is higher than observed in Musca. This is expected because the ratio decreases towards lower velocities and the 12 CO(2-1) emission in Musca also comes from non shocked gas, such that the observed ratio is a lower limit for the shock models. In red, we show the integrated brightness of 12 CO(4-3) for the same shock models, together with the 12 CO(4-3) integrated brightness interval of the blueshifted component in Musca. The observed values in Musca show a relatively good agreement with the predicted 12 CO(4-3) integrated brightness for a pre-shock density of n H 2,0 ∼ 5·10 2 cm −3 . a value around 25 • is put forward for the angle between the magnetic field and the POS. This results in a shock velocity around 0.9 km s −1 . This fits nicely with the proposed presence of slowtype low-velocity shocks, which need to have velocities between the sound speed (c s ) and v A cos(θ). With v A the alfvén speed and θ the angle between the magnetic field and the propagation direction of the shock . The expected alfvén speed towards the ambient cloud in Musca is of the order of 2.2 km s −1 , assuming a 33µG magnetic field strength at n H 2 ∼ 4·10 2 cm −3 (Crutcher 2012). We thus find that the estimated shock velocity fits with the conditions that allow a slowtype shock. An angle < 10 • between the magnetic field and the plane of the sky is required for a shock velocity > 2.2 km s −1 .
Summarising, the CO brightness temperature ratios that we observe in Musca are consistent with the predictions from low-velocity shock models (e.g. Pon et al. 2012;. Low-to mid-J CO lines may thus indeed be used to detect low-velocity shocks but more observational studies are required to better understand the quantitative contribution from shock excitation. Since SiO transitions can be used as a tracer of certain shocks (e.g. Schilke et al. 1997;Gusdorf et al. 2008a,b), it is worth noting the SiO(5-4) line was not detected towards the Musca filament. With RADEX this allows to put a 3σ upper limit on the SiO column density of ∼ 10 14 cm −2 at the filament crest, assuming n H 2 = 10 4 cm −3 . This implies that less than 0.06 % of Si is in the form of SiO at the Musca filament crest (see App. C for more information). Taking that all the gas-phase silicon is found in the form of SiO at the crest (Louvet et al. 2016), suggests that all Si is locked in the grains of the Musca cloud at an early stage of evolution. To explain the presence of SiO in the gas phase, as is found in regions like W43 (Louvet et al. 2016), implies that there is the need for an important dynamical or radiative history in such clouds to release Si from the grains in the gas phase.

A filament accretion shock signature in simulations
From the analysis of the observations, we found the presence a of warm gas layer around the dense gas (filament crest and strands) in the Musca cloud. These findings are here compared with synthetic observations from a converging flow simulation of filament formation (Clarke et al. 2018). We emphasize that this is a qualitative comparison to investigate the role of lowvelocity shocks in filament formation, because the simulation was set up for simulating slightly denser features such as the filaments in Taurus. Nevertheless, we find in the simulations some generic features related to filament accretion that are observed in Musca. In these simulations a radially convergent flow forms dense filaments with n H 2 10 4 cm −3 , see Fig. 11 and Fig. 13. The filamentary structure, that continuously accretes mass from this converging flow, has a size of 3-4 pc. Turbulence in the flow leads to structure formation in the converging flow, causing inhomogeneities in the accretion and subsequent substructure in the filament. The simulation includes hydrodynamics, self-gravity, and heating and cooling coupled with non-equilibrium chemistry. From this self-consistent CO formation in the simulations, synthetic observations are produced. The filament edge is defined by filament accretion shocks, see for example Fig. 11, which is a low-velocity shock due to mass accretion on an interstellar filament from a converging flow in a molecular cloud. For more information on these simulations and synthetic observations, we refer to Clarke et al. (2017Clarke et al. ( , 2018. No magnetic field is included in the simulations, which could have an impact on the shock properties (e.g. Draine et al. 1983;Lesaffre et al. 2013). For a shock propagating perpendicular to the magnetic field, this makes the shock wider and decreases the peak temperature. However, slow-type shocks, which can fit the observed emission, tend to resemble non-magnetised shocks because the dynamics is driven by gas pressure . In the simulations the warm post-shock layer is not resolved, as post-shock cooling of a low-velocity J-type shock occurs over small distances (e.g. Lesaffre et al. 2013;Whitworth & Jaffa 2018), implying that the gas heated by the shocks is smeared out. Consequently, the temperature of the heated gas is underestimated. Nonetheless, it is generally observed in the simulations that the 3D dense filamentary structure is confined by a gas layer that shows an increase in temperature as a result of filament accretion shocks. Figure 12 shows a temperature map, cut out of the simulation, that illustrates this temperature increase at the border of the filament, followed by a temperature drop into the filament.
The good correspondence with synthetic spectra derived from the simulation is remarkable, see Fig. 11. We observe the bulk emission of the filament with the C 18 O(2-1) line and, more importantly, a bright blueshifted component in 12 CO(4-3) with no C 18 O emission that connects the velocity of the inflowing cloud and the dense gas. This is similar to what is found in the Musca filament. In the simulations, this 12 CO(4-3) component corresponds to a strong velocity transition and density increase from a filament accretion shock of the inflowing mass reservoir, see Fig. 11. This accretion shock is responsible for the observed blueshifted 12 CO(4-3) component. It can be wondered how the slightly lower temperatures in the simulation compare to the estimated values for Musca give rise to the same signature. But as mentioned before, the simulations were designed for somewhat higher densities than Musca, which increases the 12 CO(4-3) emission and thus compensates for the slightly lower temperatures.
Since the dense filaments in the simulation are confined by an accretion shock, it is not surprising that the 12 CO(4-3) excess is found at both observed locations in Musca as it likely confines the entire filament. Figure 12 displays that behind the heated gas layer by the filament accretion shock there is also a sharp decrease in temperature from ∼15 K to ∼10 K. This creation of cold and dense gas through a filament accretion shock can also be observed in the temperature histograms of the simulations where the dense gas, created by an accretion shock, is either warmer (>25 K) or significantly colder (∼ 10 K) than the inflowing cloud (∼ 15-17 K), see Fig. 13. This demonstrates that accretion shocks from a continuously inflowing cloud, as observed in Musca, play an essential role in the cooling of the dense star forming ISM. This cooling is related to the increased shielding of the dense gas and the increased cooling rate by CO due to the higher density, since CO is an important coolant of molecular clouds (e.g. Goldsmith & Langer 1978;Whitworth & Jaffa 2018).

Filament formation due to the dissipation of MHD flows
We have put forward that the 12 CO(4-3)/ 12 CO(2-1) and 12 CO(4-3)/ 12 CO(3-2) excess in the blueshifted component can be explained by the presence of warm, dense gas that has properties predicted by low-velocity shock models. In particular the observed emission can fit with slow-type shocks, while simulations show that these slow-type shocks typically occur at relatively low velocities of v s < 3 km s −1 found in Musca Park & Ryu 2019). Furthermore, as the flow near the filament is expected to be well aligned with the magnetic field (Cox et al. 2016), the shock velocity upper limit for slow-type shocks (v A cos θ) reaches a maximum. This increases the possibility for slow-type shocks to occur. A small, but significant, change in orientation of the magnetic field takes place from the mass reservoir to the filament in Musca (Planck Collaboration et al. 2016) which could be the result of low-velocity shocks (Fogerty et al. 2017). Lastly, the observed warm and dense gas layer connects the velocity of the ambient cloud and dense gas in Musca, which indicates that the accretion shock dissipates the kinetic energy of the supersonic converging flow and increases the cooling rate to form a (trans-)sonic dense and cold filament. If the Musca filament is indeed confined by the accretion shocks, we can do a first check whether the current results give a reasonable estimated mass accretion rate. The mass accretion rate through this low-velocity shock is estimated in Appendix D, giving a mass accretion rate 20 M pc −1 Myr −1 . This is similar to the typical values reported for mass inflow towards filaments that will form low-mass stars (e.g. Palmeirim et al. 2013), and would imply that Musca can accrete its current mass in 1 Myr.

Summary
In this paper we present observations of [C II] and [O I] with SOFIA, and observations of 12 CO(4-3), 12 CO(3-2) and [C I] with APEX towards the Musca filament. Studying these important cooling lines, we find that Musca has an extremely weak FUV field (< 1 G 0 ), confirmed by determining the FUV field with a census of ionising stars. We further estimated a density n H 2 ∼  roughly indicating the high column density filament. Note the local increase in temperature at the sub-filament border due to the accretion shocks and the temperature drop in the sub-filaments after the shock compared to the inflowing mass reservoir. Because of smoothing, the temperature peaks of the shocked gas in the map have decreased. The area of the simulation that is covered by this figure is indicated in Fig.  11 5·10 2 cm −3 for the ambient gas near the Musca filament, and found that 20 to 50% of carbon in the ambient cloud is still in atomic form. In the interface region between the ambient cloud and dense gas in the crest and strands, blueshifted excess emission was found with the 12 CO(4-3), (3-2), and (2-1) lines. A non-LTE radiative transfer study with RADEX points to the existence of a small warm ( 50 K) gas layer at relatively high density (n H 2 ≥ 3·10 3 cm −3 ). This layer can actually be clumpy with pockets of warm gas embedded in a cooler interclump phase. This excess emission fits with predictions for nonirradiated slow-type low-velocity shock models with a pre-shock density n H 2,0 ∼ 5·10 2 cm −3 . The estimated shock velocity also fits with the shock velocity interval for slow-type shocks.
Comparing the location of the shock emission in the spectra with synthetic observations in simulations of filament formation, we obtain a scenario where the dense Musca filament is formed by low-velocity filament accretion shocks in a colliding flow. These low-velocity MHD filament accretion shocks dissipate the supersonic kinetic energy of the colliding flow, creating a (trans-)sonic dense filament. Because of increased shielding and efficient CO cooling in the dense post-shock gas, these accretion shocks play an important role in cooling the dense filamentary ISM that can form stars in the near future.
A&A proofs: manuscript no. template Fig. 13. Left: Density distribution of the gas with a sufficient CO abundance to be observable in a 12 CO(4-3) excess component. CO is considered detectable if more than half of the carbon is locked in CO. Note that the density distribution is relatively bimodal, with the densest already shocked gas highlighted in red. Right: Temperature distribution of all the CO gas in grey, with in red the temperature of the high density gas highlighted on the left. We find for the dense gas that almost all gas at T = 15 -25 K has disappeared, leaving only the cold (< 15 K) and warm gas (> 25 K). This demonstrates that post-shock gas is either warm or cooled behind the shocked gas layer compared to the temperature of the inflowing mass reservoir. ation via the exchange reaction between 13 C + and 12 CO leads to an enhancement of the 13 CO abundance (e.g. Liszt & Ziurys 2012;Röllig & Ossenkopf 2013;Szűcs et al. 2014). In addition -though probably less important for Musca because of the low FUV field -the more abundant 12 CO and 13 CO isotopologues shield themselves from the destructive effect of FUV photons more efficiently than the less abundant C 18 O isotopologue because the photodissociation of CO is governed by line absorption. In this context we note that Hacar et al. (2016) found indications of strong fractionation leading to a significant increase of the 13 CO/C 18 O abundance ratio above the standard value of 7 at A V < 3. Further analysis of this abundance ratio in our compan-ion paper on the Musca filament confirms this result presented in Hacar et al. (2016).
In the previous two model fits, we determined the background and foreground components which we now scale up by the abundance ratio 12 CO/ 13 CO ∼ 60. Note that this value can be significantly lower because of fractionation (see above). The resulting fit is shown in the left lower panel of Fig. B.1. Due to the large optical depths, we observe flat-topped spectra in the background and foreground layer. Since the observed intensities do not fit the natural carbon abundance ratio for the second component we only fix the first component and allowed the model to vary over a range of values in the second component to find the best possible solution. In addition we included a foreground component, which would only be observable in the 12 CO isotope, to remove the emission excess in the blueshifted wing and possibly reproduce the bump located in the blueshifted part of the observed 12 CO spectrum. The resulting spectra show a flat spectrum at velocities with bright C 18 O(2-1) emission, followed by an emission peak at the blue-shifted wing. Thus, the model results do not manage to represent the line shape, the observed bump in particular, and intensity by taking self-or foreground absorption effects of the scaled up background components derived from the isotopes, see lower left panel of Fig. B.1.
Since two background components do not represent the observed spectrum, we include an additional component to the model, which is only visible in the 12 CO spectrum but vanishes in the noise for the weaker isotopes, see lower right panel of Fig.  B.1. To meet this criteria we need to include a background component with low optical depth, which results in an increased temperature since a low optical depth requires an increase in temperature to produce the same amount of emission. This third background component is located at v = 2.8 km s −1 , see Table B.1. In addition we allow cold optically thin components in order to correct for the emission excess in the central and red-shifted part of the spectrum. As in the previous model fit we do not fix the central second component by the carbon abundance ratio, but allow it to vary in a wide range to find the best model fit to the observed data. We thus manage to represent the spectrum well, only with three background components.   environments with a radiative or dynamic history like W43 or protostellar outflows to release SiO in the gas phase.