D Arzoumanian
email: doris.arzoumanian@ias.u-psud.fr

Ph André
email: pandre@cea.fr

N Peretto

V Könyves

Formation and evolution of interstellar filaments Hints from velocity dispersion measurements ⋆,⋆⋆

Keywords: formation -ISM, clouds -ISM, structure -evolution -submillimeter, ISM

We investigate the gas velocity dispersions of a sample of filaments recently detected as part of the Herschel Gould Belt Survey in the IC 5146, Aquila, and Polaris interstellar clouds. To measure these velocity dispersions, we use 13 CO, C 18 O, and N 2 H + line observations obtained with the IRAM 30 m telescope. Correlating our velocity dispersion measurements with the filament column densities derived from Herschel data, we show that interstellar filaments can be divided into two regimes: thermally subcritical filaments, which have transonic velocity dispersions (c s < ∼ σ tot < 2 c s ) independent of column density and are gravitationally unbound; and thermally supercritical filaments, which have higher velocity dispersions scaling roughly as the square root of column density (σ tot ∝ Σ 0 0.5 )and which are self-gravitating. The higher velocity dispersions of supercritical filaments may not directly arise from supersonic interstellar turbulence but may be driven by gravitational contraction/accretion. Based on our observational results, we propose an evolutionary scenario whereby supercritical filaments undergo gravitational contraction and increase in mass per unit length through accretion of background material, while remaining in rough virial balance. We further suggest that this accretion process allows supercritical filaments to keep their approximately constant inner widths (∼0.1 pc) while contracting.

Introduction

Interstellar filaments have recently received special attention, thanks to the high quality and dynamic range of Herschel 1 observations (André et al. 2010;Men'shchikov et al. 2010;Molinari et al. 2010;Henning et al. 2010;Hill et al. 2011). The submillimeter dust continuum images of nearby molecular clouds (MCs) taken with the SPIRE (Griffin et al. 2010)a n dP A C S (Poglitsch et al. 2010) cameras onboard Herschel provide key information on both dense cores on small scales (<0.1 pc) and the filamentary structure of the parent clouds on large scales (>1 pc), making it possible to investigate the physical connection between these two components of the cold interstellar medium (ISM) and to draw a more complete picture of star formation.

Already before Herschel, filamentary structures were known to be present in MCs (e.g., Schneider & Elmegreen 1979;Myers 2009) and appeared to be easily produced by any numerical simulation of MC evolution that includes hydrodynamic (HD) or magnetohydrodynamic (MHD) turbulence (e.g., Federrath et al. 2010;Hennebelle et al. 2008;Mac Low & Klessen 2004). Individual interstellar filaments have also been recently studied with molecular line observations (e.g., Pineda et al. 2011;Hacar & Tafalla 2011;Miettinen 2012;Bourke et al. 2012;Li & Goldsmith 2012). Herschel observations now demonstrate that ⋆ Based on observations carried out with the IRAM 30 m telescope. IRAM is supported by INSU/CNRS (France), MPG (Germany), and IGN (Spain). ⋆⋆ Appendix A available in electronic form at http://www.aanda.org 1 Herschel is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA (Pilbratt et al. 2010).

these filaments are truly ubiquitous in MCs and have a direct link to the formation process of prestellar cores (André et al. 2010;Palmeirim et al. 2013).

Filaments are detected with Herschel both in star-forming regions (Könyves et al. 2010;Bontemps et al. 2010), where they are associated with the presence of prestellar cores and protostars (André et al. 2010), and in non-star-forming regions, such as the Polaris translucent cloud (Men'shchikov et al. 2010;Miville-Deschênes et al. 2010;Ward-Thompson et al. 2010). These findings support the view that filament formation precedes any star-forming activity in MCs.

In a previous work [START_REF] Arzoumanian | [END_REF], we characterized the physical properties of the filaments detected in three regions observed as part of the Herschel Gould Belt survey (HGBS, André et al. 2010) and found that all filaments share a roughly constant inner width of ∼0.1 pc regardless of their central column density and environment. This result has been confirmed in other fields of the HGBS (Peretto et al. 2012;Palmeirim et al. 2013).

Observationally, only thermally supercritical filaments for which the mass per unit length exceeds the critical value (M line > M line,crit ) show evidence of prestellar cores and the presence of star formation activity, whereas thermally subcritical filaments (M line < M line,crit ) are generally devoid of Herschel prestellar cores and protostars (André et al. 2010). The critical mass per unit length, M line,crit = 2c 2 s /G (where c s is the isothermal sound speed, and G is the gravitational constant), is the critical value required for a filament to be gravitationally unstable to radial contraction and fragmentation along its length (Inutsuka & Miyama 1997). Remarkably, this critical line mass M line,crit only depends on gas temperature (Ostriker 1964). It is approximately equal to ∼16 M ⊙ /pc for gas filaments at T = 10 K, corresponding Article published by EDP Sciences A119, page 1 of 15 to c s ∼ 0.2k ms -1 . Assuming that interstellar filaments have Gaussian radial column density profiles, an estimate of the mass per unit length is given by M line ≈ Σ 0 × W fil where W fil is the typical filament width (see Appendix A of André et al. 2010 and[START_REF] Arzoumanian | [END_REF]) and Σ 0 = µm H N 0 H 2 is the central gas surface density of the filament.

For a typical filament width of ∼0.1 pc, the theoretical value of M line,crit ∼ 16 M ⊙ /pc (for T ≈ 10 K) corresponds to a central column density N 0 H 2 ∼ 8 × 10 21 cm -2 (or visual extinction A 0

V ∼ 8). To first order, therefore, interstellar filaments may be divided into thermally supercritical and thermally subcritical filaments depending on whether their central column density N 0 H 2 is higher or lower than ∼8 × 10 21 cm -2 , respectively. Thermally supercritical filaments, with average column densities > ∼ 8 × 10 21 cm -2 , are expected to be globally unstable to radial gravitational collapse and fragmentation into prestellar cores along their lengths. The prestellar cores formed in this way are themselves expected to collapse locally into protostars. Interestingly, the critical column density ∼8 × 10 21 cm -2 of interstellar filaments is very close to the star formation threshold at a background extinction A V ∼ 6-9a b o v ew h i c ht h e bulk of prestellar cores and young stellar objects are observed in nearby clouds (e.g., Johnstone et al. 2004;Goldsmith et al. 2008;Heiderman et al. 2010;Lada et al. 2010;André et al. 2010[START_REF] André | Computational Star Formation[END_REF].

The results summarized above on the properties of interstellar filaments were made possible thanks to the sensitivity and resolution of the HGBS observations, which allowed us to detect structures down to A V ∼ 0.1 and to resolve individual filaments in nearby regions (d < 0.5 kpc), down to a typical Jeans length of ∼0.02 pc. While the Herschel images already are a powerful tool to assess the importance of different physical processes involved in the formation and evolution of interstellar filaments, they do not provide any information on the underlying kinematics of the gas.

In the present work, we complement our first results obtained from Herschel imaging data with molecular line observations to gain insight into the velocity dispersion of the interstellar filaments and reinforce our understanding of filament formation and evolution. We present 13 CO, C 18 O, and N 2 H + line observations, taken with the IRAM 30 m-telescope, toward a sample of filaments detected with Herschel in the IC 5146, Aquila Rift, and Polaris flare regions. In Sect. 2 we present the observations and the data reduction. Sections 3 and 4 summarize the data analysis, Sect. 5 presents the results and Sect. 6 discusses their implications in the proposed scenario for interstellar filament formation and evolution.

IRAM observations and data reduction

Molecular line observations of a sample of filaments in IC 5146, Aquila, and Polaris were carried on at the IRAM-30 m telescope at Pico Veleta, Spain, during two observing runs in March and August 2011.

In March 2011, we observed the C 18 O(2-1) and 13 CO(2-1) spectral lines with the HEterodyne Receiver Array (HERA, Schuster et al. 2004)a n dt h eN 2 H + (1-0) line with the Eight MIxer Receiver (EMIR, a multiband mm-wave receiver, Carter et al. 2012) toward the filaments of IC 5146. In addition to single-position spectra taken in the on-off observing mode, we mapped filament 14 in IC 5146 (yellow rectangle in Fig. 1)inthe C 18 O(2-1) and 13 CO(2-1) transitions using the on-the-fly mapping mode with HERA. The Aquila rift and the Polaris clouds The spatial resolutions of the observations correspond to the 12 ′′ ,2 3 ′′ and 28 ′′ FWHM beamwidths of the 30 m antenna at 219.6 GHz [C 18 O(2-1)], 109.8 GHz [C 18 O(1-0)], and 93.2 GHz [N 2 H + (1-0)], respectively. We used the VESPA autocorrelator as backend with a spectral resolution of 20 KHz. This spectral resolution translates to velocity resolutions of 0.027, 0.055, and 0.064 km s -1 at the frequencies of the C 18 O(2-1), C 18 O(1-0), and N 2 H + (1-0) transitions, respectively.

We used the frequency-switching mode for the CO observations and the position-switching mode for the N 2 H + observations, using a reference position offseted by ∼10 ′ -15 ′ from each target position. The off position was selected from the Herschel images and the lack of N 2 H + emission was checked using observations in the frequency-switching mode. During the observations, calibration was achieved by measuring the emission from the sky, an ambient load and a cold load every ∼15 min. The telescope pointing was checked and adjusted every ∼2h and the pointing accuracy was found to be better than 3 ′′ .T h e typical rms of each single position (on-off mode) spectrum is ∼0.05 K (in T * A scale). All of the data were reduced with the GILDAS/CLASS software package (the Grenoble Image and Line Data Analysis System, a software provided and actively developed by IRAM 2 .

Filament sample

A total number of 100 positions toward 70 filaments or segments of filaments3 were observed in the three regions (IC 5146,Aquila,Polaris). Most of the densest filaments, with N H 2 > ∼ 1 × 10 22 cm -2 , were detected in N 2 H + (1-0), which is a low Fig. 2. Spectra observed in the C 18 O(2-1) and the N 2 H + (1-0) transitions, toward filament 14 (top)a n d2( bottom), respectively. The best Gaussian fits to the spectra are overlaid in blue. The fitted parameters are given in Table 1. optical depth tracer of dense gas, as well as in C 18 Oa nd 13 CO. The lower column density filaments were only detected in C 18 O and 13 CO, while the faintest filaments were not detected 4 .

The final sample discussed in the remainder of this paper corresponds to the 46 filaments detected in either N 2 H + or C 18 O. This sample consists of 44 individual filaments detected with the IRAM 30 m telescope in IC 5146 (11 filaments), Aquila rift (26 filaments), and Polaris flare (7 filaments), plus 2 filaments corresponding to the NGC 2264C elongated clump (Peretto et al. 2006) and the DR21 ridge (Schneider et al. 2010;Hennemann et al. 2012).

The positions of the filaments observed in IC 5146 are shown in Fig. 1, overlaid on the Herschel column density map. The filaments observed in Aquila and Polaris are shown in Figs. A.5, A.8 and A.10, respectively. A N 2 H + (1-0) spectrum and a C 18 O(2-1) spectrum toward two filaments observed in IC 5146 are shown in Fig. 2. The positions and derived parameters of all the filaments used for the analysis presented in this paper are summarized in Table 1 and the corresponding spectra are shown in online Appendix A.

Line-of-sight velocity dispersions

In order to combine the velocity information obtained with various line tracers, we estimated the total velocity dispersion of the mean free particle in molecular clouds. 4 A few filaments observed in IC 5146 (5 filaments) and Polaris (1 filament) were only detected in 13 CO. The velocity dispersions (∼0.3 km s -1 on average) of these filaments were not included in the filament sample due to larger uncertainties in the optical depth of this transition.

We estimated the non-thermal velocity dispersion of the gas by subtracting the thermal velocity dispersion from the linewidth measured for each species, assuming that the two contributions are independent of each other and could be added in quadrature (Myers 1983). The thermal velocity dispersion of each species is given by

σ T (µ obs ) = k B T µ obs m H , (1)
where k B is the Boltzmann constant, µ obs is the atomic weight of the observed molecule, i.e. µ obs = 29 for N 2 H + and µ obs = 30 for C 18 O. We adopted a gas temperature of 10 K which is the typical temperature of starless molecular clouds and is close to the dust temperature measured with Herschel for the present sample of filaments. The non-thermal velocity dispersion is then equal to

σ NT = σ 2 obs -σ 2 T (µ obs ), (2)
where σ obs =∆ V/ √ 8ln2 and ∆V is the measured FWHM linewidth of the observed spectra. The linewidths were estimated by single Gaussian fits to the C 18 O spectra (cf. Appendix A for more details) and by multiple Gaussian fits to the seven components of the hyperfine multiplet of the N 2 H + transition, using the Gaussian HyperFine Structure (HFS) fitting routine of the CLASS software package. This routine derives the line optical depth by assuming the same excitation temperature for all hyperfine components, and therefore yields an estimate of the intrinsic linewidth.

The total velocity dispersion of the mean free particle of molecular weight µ = 2.33 is finally given by:

σ tot (µ) = σ 2 NT + σ 2 T (µ), (3)
where

σ T (µ) = k B T
µm H ∼ 0.2k ms -1 for T = 10 K. The typical gas temperature in the center of the dense filaments of our sample is approximately T gas ∼ 10 K since the dust temperature derived from Herschel data is T dust ∼ 10 K and T gas is expected to be well coupled to T dust at densities > ∼ 3 × 10 4 cm -3 (see Galli et al. 2002). As to the lower-density filaments in our sample, CO observations directly show that their gas temperature is also close to 10 K (e.g., Goldsmith et al. 2008;Heyer et al. 2009).

Velocity dispersion along filament 14 in IC 5146

The velocity dispersions of the filaments in our sample have been estimated from single-position observations toward each filament. To assess the reliability of inferring the velocity dispersion of a filament from a single position spectrum, we investigated the variations of the velocity dispersion along the crest of filament 14 in IC 5146 which we mapped in both 13 CO(2-1) and C 18 O(2-1) -see Fig. 3.

From the 13 CO(2-1)/C 18 O(2-1) ratio map, we first estimated the optical depth to check whether or not the C 18 O(2-1) emission is optically thin.

Assuming uniform excitation temperature along the line of sight, 13 CO(2-1) and C 18 O(2-1) line emissions can be expressed as a function of optical depth as follows:

T * ,i A (υ) = T i ex [1 -exp(-τ i (υ))], (4)
where i corresponds to either C 18 O(2-1) or 13 CO(2-1), T * ,i A (υ) is the detected antenna temperature for each velocity channel υ, A119, page 3 of 15 A&A 553, A119 (2013) Fig. 3. Left: Herschel column density map of filament 14 in IC 5146 corresponding to the yellow rectangle in Fig. 1, in units of visual extinction (A V )whereN H 2 ∼ A V × 10 21 cm -2 . The resolution of the map is 36.9 ′′ . The black contours correspond to column densities of 1 and 2 × 10 21 cm -2 . These contours are reproduced in the middle and right panels. Middle: 13 CO(2-1) integrated intensity map (in units of K km s -1 ) over the LSR velocity range -2 .5to1k ms -1 . The map has been smoothed to a resolution of 18.5 ′′ , which corresponds to half the resolution of the Herschel column density map. The 13 CO(2-1) integrated emission of the filament is concentrated within the 10 21 cm -2 column density contour, which traces the elongated structure of the filament. Right:sameasthe middle panel for the C 18 O(2-1) transition.

T i

ex is the excitation temperature which we assume to be the same for the two isotopomers and τ i 21 is the optical depth of the (2-1) transition of C 18 Oor 13 CO.

The optical depth values at the systemic velocity (υ sys )o f the cloud for each transition were estimated from the ratio of the 13 CO(2-1) and C 18 O(2-1) line intensities:

T * ,13 A (υ sys ) T * ,18 A (υ sys ) = T 13 ex [1 -exp(-τ 13 21 (υ sys ))] T 18 ex [1 -exp(-τ 18 21 (υ sys ))] , (5)
where T 13 ex = T 18 ex and τ 13 21 = Xτ 18 21 , with X = [ 13 CO]/[C 18 O] = 5.5 (the assumed mean value of the abundance ratio in the local ISM, e.g. Wilson & Rood 1994).

Maps of the mean optical depth of filament 14 were constructed from the observed integrated intensity ratio map of the 13 CO(2-1) and C 18 O(2-1) transitions, over the velocity range -2.5 to +1kms -1 ,

R 13/18 = 1 -2.5 T * ,13 A (υ)dυ/ 1 -2.5 T * ,18 A (υ)dυ.
We estimated the mean values of τ 18 21 and τ 13 21 (corresponding to the optical depth values integrated over the velocity range -2.5 to +1k m s -1 ) from the R 13/18 ratio map by solving a velocity-integrated form of Eq. ( 5)o nap i x e lb yp i x e lb a s i s . Figure 4 shows the mean optical depth maps of the 13 CO(2-1) and C 18 O(2-1) transitions for filament 14 (middle and right panels, respectively) along with the R 13/18 ratio map (left panel).

While the bulk of the C 18 O emission is optically thin (τ 18 21 < 0.6) along the main axis of the filament, the 13 CO emission is significantly optically thick (τ 13 21 > 2) in the southern half part of the filament and only marginally optically thin (τ 13 21 ∼ 0.5-2) in the northern part of the filament.

Figure 5 shows the positions along filament 14 where spectra were extracted from the maps (left-hand side of the figure) and the velocity dispersions measured along the filament (right). The velocity dispersions derived from the 13 CO(2-1) spectra were corrected for the broadening effect due to the finite optical depth of the 13 CO(2-1) transition. To do so, we measured the variations of the FWHM linewidth, ∆V obs , of a simple Gaussian line model, T model τ 0 (υ), as a function of optical depth. The model line spectrum had the form T model

τ 0 (υ) = T ex [1 -exp -τ(υ) ], where T ex is a uniform excitation temperature, τ(υ) = τ 0 exp -4ln2(υ-υ 0 ) 2
∆V 2 int corresponds to an intrinsically Gaussian distribution of optical depths with velocity, υ,a n d∆V int is the intrinsic FWHM linewidth of the model. By fitting Gaussian profiles to synthetic spectra, we estimated the FWHM width ∆V obs of the model as a function of τ 0 , and derived correction factors ∆V int /∆V obs ≈ 0.77-0.92, for 13 CO(2-1) optical depths ranging between 0.5 and 2.5. The most extreme correction would be a factor of 0.74 for a 13 CO(2-1) optical depth of 3, if the [ 13 CO]/[C 18 O] abundance ratio has a value of 7 (e.g., Wu et al. 2012), i.e., higher than the standard value in the local ISM.

Inspection of the linewdiths measured along filament 14 shows that the variations of the internal velocity dispersion are A (υ)dυ/ 1 -2.5 T * ,18 A (υ)dυ for filament 14 in IC 5146. Only those pixels where the column density derived from Herschel is higher than10 21 H 2 cm -2 have been considered (black contour in Fig. 3). Middle: map of the mean optical depth of the 13 CO(2-1) transition. The contours in green correspond to a value of 1. Right: map of the mean optical depth of the C 18 O(2-1) transition. The mean optical depth of the C 18 O transition is τ(C 18 O) < 1 over the entire filament. 6). Right: total velocity dispersions, σ tot (see Eq. ( 3)), measured along thesamefilament.Thex-axis goes from the south to the north of the filament. The red triangle and blue squares correspond to measurements made on the C 18 O(2-1) and 13 CO(2-1) maps, respectively. Individual spectra were constructed by spatially averaging the signal located within 10 ′′ of the central position to increase the signal to noise ratio of the data. The linewidths of the 13 CO(2-1) spectra were corrected for the broadening due to optical depth effects (see explanations in Sect. 4). The median value of σ tot is 0.25 km s -1 with a narrow dispersion of 0.03 km s -1 .

small along the crest of the filament (cf. right panel of Fig. 5), with an average velocity dispersion (0.25 ± 0.03) km s -1 .This is consistent with previous observations of a few low-density filaments showing that these structures have velocity dispersions close to the thermal velocity dispersion ∼0.2 km s -1 for T = 10 K (cf. Hily-Blant 2004; Hily-Blant & Falgarone 2009) and which do not vary much along their length (Hacar & Tafalla 2011;Pineda et al. 2011). Recently, Li & Goldsmith (2012) studied the Taurus B213 filament and found that it is characterized by a coherent velocity dispersion of about ∼0.3 km s -1 . These results suggest that the velocity dispersion observed at a single position toward a filament provides a reasonably good estimate of the velocity dispersion of the entire filament 5 . Nevertheless mapping observations of a broader sample 5 An interstellar filament is an elongated structure characterized by small column density variations along its crest (less than a factor of 3; see Fig. 3-left). Therefore, the velocity dispersion variations induced by the trend σ tot ∝ Σ 0 0.5 found in Sect. 5 below for supercritical filaments remain small (≪2) along a given filament.

A119, page 5 of 15 A&A 553, A119 (2013) of filaments in several regions would be valuable to confirm the robustness of this assumption.

(J2000) (J2000) [10 21 cm -2 ] M ⊙ /pc [km s -1 ][ K ( T A *)] [km s -1 ][ k m s -1 ] (1) (2) (3) (4) (5) (6) (7) (8) (

Two regimes for interstellar filaments

In this section, we investigate the variation of internal velocity dispersion with column density for our sample of filaments.

Figure 6 shows the total velocity dispersion as a function of central column density for the 46 filaments. Values for the column density and velocity dispersion of the NGC 2264C elongated clump (red square in the following plots) are taken from Peretto et al. (2006). Peretto et al. estimated the column density from 1.2 mm dust continuum mapping observations and the velocity dispersion of NGC 2264C from N 2 H + (1-0) observations, both obtained with the IRAM 30 m telescope. The velocity dispersion was averaged over the entire NCG 2264C clump and reflects the velocity dispersion of the filamentary clump as opposed to the embedded protostellar cores. The velocity dispersion of the DR21 filament (green square in the following plots) is taken from Schneider et al. (2010) who derived it from N 2 H + (1-0) observations, while the column density of DR21 comes from Herschel observations (Hennemann et al. 2012) obtained as part of the HOBYS key project (Motte et al. 2010). The black triangle corresponds to filament 14 in IC 5146 and the associated error bar comes from the dispersion of the total velocity dispersion along this filament as extracted from the map shown in Fig. 5. The vertical error bars of the other points correspond to the errors in estimating the linewidth from Gaussian fits to the observed spectra, which are on the order of ∼0.01 km s -1 . The error on the measured velocity dispersion is ∼0.03 km s -1 for the filaments observed in Aquila and IC 5146 and 0.1 km s -1 for DR21 and NGC 2264C (representing the variation of σ tot along the filament crests). The errors on the filament central column densities come from the uncertainties in deriving column densities6 from Herschel dust continuum maps and are typically a factor of ∼2 mainly due to uncertain assumptions concerning the adopted dust opacity law (see Könyves et al. 2010;[START_REF] Arzoumanian | [END_REF]. A recent study comparing Herschel and near-IR extinction data in Orion A reported a weak trend between dust opacity and column density, at least in the regime 1 < ∼ A V < ∼ 10, which has been interpreted as evidence of dust grain evolution with environment (Roy et al. 2013). But the uncertainty induced by this effect remains small for most of the filaments in our sample.

An additional systematic effect on the column density comes from the random inclination angles of the filaments to the plane of the sky, which implies an overestimation of the intrinsic column densities (and masses per unit length) by ∼60% on average for a large sample of filaments (for more details see Appendix A of [START_REF] Arzoumanian | [END_REF]. This systematic correction factor is shown as a horizontal arrow pointing to the left in Fig. 6.

The results of our velocity dispersion measurements confirm that the filaments of our sample can be divided into two physically distinct groups based on their central column density. Filaments with column densities ≤8 × 10 21 cm -2 seem to have roughly constant velocity dispersions close to or slightly larger than the sound speed (c s < ∼ σ tot < 2 c s ), while denser filaments with column densities ≥8 × 10 21 cm -2 have velocity dispersions which increase as a function of column density (cf. Fig. 6). The upper x-axis of Fig. 6 shows an approximate mass-per-unitlength scale derived from the bottom x-axis scale by multiplying the central column density of each filament by a characteristic filament width W fil ∼ 0.1 pc (i.e., M line =Σ 0 × W fil -see Sect. 1). It can be seen from Fig. 6 that the critical mass per unit length M line,crit ∼ 16 M ⊙ /pc for T ≈ 10 K introduced in Sect. 1 corresponds to a column density boundary which divides the filaments into two regimes where thermally subcritical filaments (M line < M line,crit ) have roughly constant velocity dispersions with a mean value of 0.26 ± 0.05 km s -1 , while thermally supercritical filaments (M line > M line,crit ) have velocity dispersions which increase as a function of projected column density as a Fig. 6. Total velocity dispersion versus observed central column density (uncorrected for inclination effects): blue for IC 5146, magenta for Aquila, cyan for Polaris, red for NGC 2264C, and green for DR21 filaments. The squares and triangles correspond to velocity dispersions measured from N 2 H + spectra and C 18 O spectra, respectively. See text for details on the error bars and for uncertainties in the column densities. The horizontal dashed line shows the value of the sound speed ∼0.2k ms -1 for T = 10 K. The vertical dashed line marks the boundary between thermally subcritical and thermally supercritical filaments where the estimated mass per unit length M line is approximately equal to the critical value M line,crit ∼ 16 M ⊙ /pc for T = 10 K, equivalent to a column density of 8 × 10 21 cm -2 . The grey band shows a dispersion of a factor of 3 around this nominal value. The dotted line running from the bottom left to the top right corresponds to σ tot ∝ N H 2 0.5 normalized to 0.2 km s -1 at the subcritical/supercritical boundary. The blue solid line shows the best power-law fit σ tot ∝ N H 2 0.35 ± 0.14 to the supercritical filaments.

power law σ tot ∝ N H 2 ,obs 0.35 ± 0.14 . The latter relation becomes σ tot ∝ N H 2 ,corr 0.41 ± 0.15 if "intrinsic" column density values are used (after correcting the observed values for a ∼60% overestimation on average due to inclination effects). The division of the filament sample into two regimes does not correspond to a sharp boundary but to a narrow border zone represented by the grey band in Fig. 6 which results from 1) a spread of a factor of < ∼ 2 in the effective width of the filaments around the nominal value of 0.1 pc used to convert the critical mass per unit length to a critical column density; 2) uncertainties in the intrinsic column densities of the filaments due to a random distribution of viewing angles; and 3) a small dispersion in the effective gas temperature of the filaments, implying a small spread in the theoretical value of the effective critical mass per unit length (see below).

To derive more accurate filament masses per unit length, we constructed radial column density profiles from the Herschel column density maps, taking a perpendicular cut to each filament at the position observed with the 30 m telescope. The observed mass per unit length, M obs line , of each filament was then derived by integrating the measured column density profile over radius. The mass per unit length of DR21 corresponds to the average value along the crest of the filament derived by Hennemann et al. (2012), while the value for NGC 2264C was obtained from its estimated mass (1650 M ⊙ ) and length (0.8 pc). These M obs line estimates of the mass per unit length are used in the following discussion instead of the simpler estimates shown at the top of Fig. 6.

Figure 7a shows the total velocity dispersion of each filament (same y-axis as in Fig. 6) as a function of the estimated A119, page 7 of 15 A&A 553, A119 (2013) Fig. 7. a) Total velocity dispersion versus intrinsic mass per unit length (corrected for inclination effects) for the sample of 46 filaments. The symbols and colors, as well as the y-axis of the plot, are the same as in Fig. 6. The horizontal segment at the bottom of the plot shows the typical error bar on the estimated masses per unit length. Specific error bars are shown for NGC 2264C and DR21 since they are larger 5 . The horizontal dashed line shows the value of the thermal sound speed ∼0.2kms -1 for T gas = 10 K. The vertical grey band marks the theoretical position of the effective critical mass per unit length M eff line,crit for effective temperatures in the range 10 < ∼ T eff < ∼ 20 K (see text). The dotted line running from bottom left to top right corresponds to M line = M line,vir /2, which marks the boundary between unbound and bound filaments located on the left and right hand side of the border line, respectively. The blue solid line, σ tot ∝ M 0.36 ± 0.07 line , shows the best power-law fit for the bound filaments. The best-fit exponent remains almost unchanged (0.31 ± 0.08) when the fitting is performed after removing DR21 and NGC 2264C from the sample. b) Corrected virial parameter (α line,vir ) versus intrinsic mass per unit length for the same filaments (x-axis is the same as in panel a)). The vertical grey band is the same as in panel a) and separates subcritical filaments (M obs line < M eff line,crit ) on the left from supercritical filaments (M obs line > M eff line,crit ) on the right. The horizontal dashed line (α line,vir = 2) shows the boundary between gravitationally unbound (α line,vir > 2) and bound (α line,vir < 2) filaments, where α line,vir = M line,vir /M obs line . The two blue solid lines show the best fits for the subcritical and supercritical filaments, respectively (see text for details). intrinsic mass per unit length obtained from the observed value after applying a 60% correction factor to account for the average overestimation of the intrinsic column density due to projection effects. This panel can be divided into two parts by comparing the intrinsic mass per unit length to the virial mass per unit length M line,vir = 2σ 2 tot /G (Fiege & Pudritz 2000), where σ tot is the observed total velocity dispersion (instead of the thermal sound speed used in the expression of M line,crit ). The dividing line is defined by M line,vir = 2 M obs line equivalent to α line,vir = 2, where

α line,vir = M line,vir /M obs line (6)
is the virial parameter of a filament 7 . Figure 7b shows how the virial parameter depends on the intrinsic mass per unit length for our sample of filaments. This plot (Fig. 7b) can be divided into four quadrants. The effective critical mass per unit length M eff line,crit = 2σ eff tot 2 /G divides the sample into subcritical filaments on the left-hand side and supercritical filaments on the right-hand side 8 . Conceptually, we define σ eff tot as the total velocity dispersion of a filament on the verge of global radial collapse. Observationally, the range of internal velocity dispersions σ tot measured here for thermally subcritical and nearly critical filaments suggests that σ eff tot corresponds to effective temperatures T eff ≡ µ m H σ eff tot 2 /k B in the 7 For a sphere of uniform density α vir = M vir /M where M vir ∼ 5Rσ 2 /G and M is the mass of the object (e.g., Bertoldi & McKee 1992). 8 In this paper, we use the terms "supercritical" and "subcritical" to characterize the mass per unit length of a filament, independently of the magnetic field strength (but see Fiege & Pudritz 2000 for how M line,crit is modified in the presence of magnetic fields). This should not be confused with the terms "(magnetically) supercritical and subcritical" used to characterize clouds with supercritical and subcritical mass-to-flux ratios, respectively. range 10 < ∼ T eff < ∼ 20 K. This implies an effective critical mass per unit length M eff line,crit in the range between ∼16 and 32 M ⊙ /pc, which is shown as a grey band in the plots of Figs. 7aa n db . The virial parameter, α line,vir , divides the sample of filaments into unbound (α line,vir > 2) structures in the upper part of the diagram and gravitationally bound structures (α line,vir < 2) in the lower part of the diagram. The virial parameter decreases as α line,vir ∝ M -0.95±0.12 line for the subcritical filaments, while this dependence is much shallower, α line,vir ∝ M -0.24±0.29 line , for the supercritical filaments of our sample. The best-fit power law exponent for the supercritical filaments remains almost unchanged (-0.28 ± 0.30) when the linear fitting is performed after removing DR21 and NGC 2264C from the sample. Figure 7b indicates that thermally supercritical filaments tend to be gravitationally bound with α line,vir = 1.0 ± 0.5 on average, and that, conversely, all thermally subcritical filaments are unbound with α line,vir > 2.

Discussion and conclusions

The results of our 13 CO, C 18 O, and N 2 H + line observations toward a sample of filaments previously detected with Herschel show that, to an excellent approximation, interstellar filaments may be divided into bound and unbound structures depending on whether their mass per unit length M line is larger or smaller than the thermal value of the critical mass per unit length M line,crit = 2c 2 s /G. Indeed, thermally supercritical filaments are found to be self-gravitating structures in rough virial balance with M obs line ∼ M line,vir , while thermally subcritical filaments appear to be unbound structures with M obs line ≤ M line,vir /2andtransonic line-of-sight velocity dispersions (c s < ∼ σ tot < 2 c s ). This confirms the usefulness of the simple gravitational instability criterion based on M line,crit which was adopted by André et al. (2010) in their analysis of the first results from the HGBS.

The detection, in the same cloud (e.g. Aquila), of both subcritical and supercritical filaments, with transonic and supersonic velocity dispersions respectively, suggests that the internal velocity dispersion within a filament is not directly related to the level of large-scale turbulence in the parent cloud. Instead, we propose that the internal velocity dispersion partly reflects the evolutionary state of a filament (see below).

6.1. An evolutionary scenario for interstellar filaments?

In the light of the results presented in this paper, we may revisit the issue of the formation and evolution of interstellar filaments and improve the scenario proposed in our earlier work based on HGBS observations (e.g. [START_REF] Arzoumanian | [END_REF].

The observed omnipresence of filamentary structures in gravitationally unbound complexes such as the Polaris translucent cloud (Men'shchikov et al. 2010;Miville-Deschênes et al. 2010;Ward-Thompson et al. 2010) suggests that large-scale turbulence rather than large-scale gravity play the dominant role in forming interstellar filaments. The finding of a characteristic filament width ∼0.1 pc [START_REF] Arzoumanian | [END_REF], corresponding to better than a factor of ∼2 to the sonic scale below which interstellar turbulence becomes subsonic in diffuse molecular gas (cf. Goodman et al. 1998;Falgarone et al. 2009;Federrath et al. 2010), also supports the view that filaments may form by turbulent compression of interstellar gas (Padoan et al. 2001). In such a picture, filaments coincide with stagnation gas associated with regions of locally converging turbulent motions, where compression is at a maximum and relative velocity differences at a minimum, and are thus expected to have relatively low (transonic) internal velocity dispersions (cf. Klessen et al. 2005). The line-of-sight velocity dispersions measured here for subcritical filaments, which all are in the range c s < ∼ σ tot < ∼ 2 c s , are consistent with this picture. Moreover, our 13 CO(2-1) and C 18 O(2-1) mapping results for filament 14 in IC 5146 (cf. Figs. 3 and5 see also Hacar & Tafalla 2011 for other filaments in Taurus) suggest that subcritical filaments are velocity-coherent structures in MCs, with small and roughly constant velocity dispersion along their length (cf. Goodman et al. 1998;Heyer et al. 2009,forthe definition of the concept of coherence in the context of dense cores).

Since subcritical filaments are unbound, they may be expected to widen and disperse on a turbulent crossing time. Given the mean velocity dispersion σ tot ∼ 0.3k ms -1 measured here for subcritical filaments, the typical crossing time of a ∼0.1pc-wide filament is ∼0.3 Myr. Subcritical filaments may be expected to disperse on this timescale unless they are confined by some external pressure as suggested by Fischera & Martin (2012). These authors modeled interstellar filaments as quasiequilibrium isothermal structures in pressure balance with a typical ambient ISM pressure P ext ∼ 2-5 × 10 4 Kcm -3 (Fischera & Martin 2012). The predicted widths of their model filaments are in rough agreement with the observed filament width of ∼0.1 pc. Moreover, the agreement between the expected and the observed widths improves when polytropic models of pressure-confined filaments, obeying a non-isothermal equation of state P ∝ ρ γ with γ<1 are considered (Inutsuka, priv. comm.).

The observed increase in velocity dispersion with increasing central column density for thermally supercritical filaments, roughly consistent with the scaling σ tot ∝ N H 2 0.5 (cf. Fig. 6), provides interesting clues to the evolution of dense, star-forming filaments. Since the filaments on the right-hand side of Fig. 6 have α line,vir ∼ 1 and are self-gravitating (see Sect. 5), they should be unstable to radial collapse and thus have a natural tendency tot /GΣ 0 ) versus filament intrinsic central column density. The vertical dashed line and the grey band mark the border zone between subcritical and supercritical filaments (cf. Fig. 6). The dotted line running from the top left to the bottom right corresponds to the thermal Jeans diameter (2R J ). The blue solid lines correspond to the best fits for the subcritical (D sub J,eff ∝ N H 2 -1.1 ± 0.1 ) and supercritical (D sup J,eff ∝ N H 2 -0.2 ± 0.1 ) filament subsamples, respectively. The exponent of the power-law fit for the supercritical filaments becomes -0.4 ± 0.1 after excluding DR21 and NGC 2264C from the sample. In contrast to the thermal Jeans length, the effective Jeans diameter of the supercritical filaments remains roughly constant with an average value of 0.14 ± 0.07 pc.

to contract with time (cf. Inutsuka & Miyama 1992;Kawachi & Hanawa 1998). The shape of their radial column density profiles with a high degree of symmetry about the filament major axis and a well-defined power law regime at large radii, approximately consistent with a radial density structure ρ ∝ r -2 [START_REF] Arzoumanian | [END_REF]Palmeirim et al. 2013;Hill et al. 2012), is consistent with a theoretical model of a nearly isothermal collapsing cylinder (Kawachi & Hanawa 1998;Palmeirim et al. 2013), and therefore also suggestive of gravitational contraction. If supercritical filaments are indeed contracting, then their central column density is expected to increase with time (cf. Kawachi & Hanawa 1998). It is thus tempting to suggest that the right-hand side of Fig. 6 may, at least partly, represent an evolutionary sequence for supercritical filaments, in which central column density corresponds to an evolutionary indicator, increasing with time. As the internal velocity dispersion of a filament provides a direct measure of its virial mass per unit length and supercritical filaments have M line > ∼ M line,vir , the trend seen on the right-hand side of Fig. 6 suggests that the mass per unit length of supercritical filaments increases with time through accretion of background material. More direct evidence of this accretion process and growth in mass per unit length for supercritical filaments exists in several cases in the form of low-density striations or subfilaments observed perpendicular to the main filaments and apparently feeding them from the side. Examples include the B211/B213 filament in Taurus where the typical velocities expected for the infalling material are consistent with the existing kinematical constraints from CO observations (Goldsmith et al. 2008;Palmeirim et al. 2013), the Musca filament (Cox et al. 2013, in prep.), and the DR21 ridge in Cygnus X (Schneider et al. 2010).

If such an accretion process accompanied by an increase in internal velocity dispersion does indeed occur in supercritical filaments, it may explain how these filaments can maintain a roughly constant inner width while contracting. In collapse A&A 553, A119 (2013) models of nearly-isothermal cylindrical filaments (cf. Inutsuka & Miyama 1992;Kawachi & Hanawa 1998), the flat inner portion of the radial density profile has a radius corresponding to the instantaneous central Jeans length R J ∼ c 2 s /GΣ 0 ,whereΣ 0 is again the central mass column density of the model filament. In the presence of a nonthermal component to the velocity dispersion, we may thus expect the central diameter of a contracting filament to be roughly given by twice the effective Jeans length, i.e., 2 R J,eff ≡ D J,eff ∼ 2 σ 2 tot /GΣ 0 .T h i se ffective Jeans diameter is plotted in Fig. 8 as a function of corrected central column density for the sample of filaments studied in the present paper. It can be seen in Fig. 8 that D J,eff is approximately constant for the thermally supercritical filaments in our sample, with an average value of 0.14 ± 0.07 pc.

In this scenario, the enhanced velocity dispersions measured for thermally supercritical filaments would not directly arise from large-scale interstellar turbulence but from the accretion of ambient cloud gas and the gravitational amplification of initial velocity fluctuations through the conversion of gravitational energy into kinetic energy during the accretion and contraction process (see, e.g., the numerical simulations presented by Peretto et al. 2007, to explain the observed characteristics of NGC 2264C). This would also be consistent with the concept of accretion-driven turbulence proposed by Klessen & Hennebelle (2010) to explain the origin of turbulent motions in a wide range of astrophysical objects.

At the same time as they contract and accrete material from the background cloud, supercritical filaments are expected to fragment into cores and form (proto)stars (Inutsuka & Miyama 1997;Pon et al. 2011). Indeed, for filamentary clouds, fragmentation into cores and subsequent local core collapse occur faster than global cloud collapse (Kawachi & Hanawa 1998;Pon et al. 2011;Toalá et al. 2012). As pointed out by Larson (2005), the filamentary geometry is thus especially favorable for fragmentation and core formation. This is also supported by observations of supercritical filaments which are usually found to harbor several prestellar cores and Class 0/Class I protostars along their length (cf. André et al. 2010). Core formation is particularly well documented in the case of the B211/B213 filament in Taurus (Schmalzl et al. 2010;Onishi et al. 2002)ortheSerpens South filament in the Aquila Rift (Gutermuth et al. 2008;Maury et al. 2011).

To assess the reliability of this tentative scenario for filament formation and evolution, comparison with dedicated numerical simulations would be very valuable. More extensive molecular line mapping of a larger sample of filaments are also very desirable, in order to set stronger observational constraints on the dynamics of these structures. André et al. (2010). The curvelet component is obtained using a morphological component analysis algorithm (MCA, from Starck et al. 2003) which enhances the contrast of the filamentary structure against the more diffuse background of the cloud. The positions of the observed spectra are plotted in red squares and blue triangles for N 2 H + and C 18 O, respectively. The numbers correspond to the filaments listed in Table 1 and are given at the upper right of the spectra in 2010). The positions of the observed spectra are plotted in red squares and blue triangles for N 2 H + and C 18 O, respectively. The numbers correspond to the filaments listed in Table 1. The corresponding Gaussian fits are highlighted in blue. Note 1: the spectrum observed toward filament 4 shows two velocity components. Both velocity components have similar linewidths (0.19 km s -1 and 0.14 km s -1 , respectively), they are separated by ∼0.8 km s -1 . They most probably belong to two filaments which are seen on the same line-of-sight (but a large map would be needed to study the kinematics of the field). The component observed at -4.5 km s -1 has been taken to be representative of filament 4.

Fig. 1 .

 1 Fig. 1. Herschel column density map of the northern and southern streamers of IC 5146 with the network of filaments highlighted in blue (adapted from Fig. 3b of Arzoumanian et al. 2011). The positions of the observed spectra are plotted as red squares and blue triangles for N 2 H + and C 18 O, respectively (the corresponding spectra of the positions observed for filament 2 and filament 14 are shown in Fig. 2). The yellow rectangle shows the subcritical filament 14 which was mapped in C 18 O(2-1) and 13 CO(2-1); cf. Fig. 3.

Fig. 4 .

 4 Fig. 4. Left: integrated intensity ratio map 1 -2.5 T * ,13A (υ)dυ/

Fig. 5 .

 5 Fig. 5. Left:integrated 13 CO(2-1) and C 18 O(2-1) intensity maps of filament 14 in IC 5146 (cf. Fig. 3). Two 13 CO(2-1) contours are overlaid: 1.2 and 1.5 K km s -1 . The blue squares and red triangles correspond to the positions along the filament where the line-of-sight velocity dispersion has been measured from the 13 CO(2-1) and C 18 O(2-1) maps respectively. The two black triangles on the C 18 O(2-1) map correspond to the positions of the central pixel of the array receiver of HERA observed in a single pointing. The labels indicate the offsets in arcsec relative to the center of the map at Ra(J2000) = 21:47:56 and Dec(J2000) = 48:08:55. The black triangle at position (-94 ′′ , -223 ′′ ) is the point which represents filament 14 in the column density vs. velocity dispersion diagram (black triangle in Fig.6). Right: total velocity dispersions, σ tot (see Eq. (3)), measured along thesamefilament.Thex-axis goes from the south to the north of the filament. The red triangle and blue squares correspond to measurements made on the C 18 O(2-1) and 13 CO(2-1) maps, respectively. Individual spectra were constructed by spatially averaging the signal located within 10 ′′ of the central position to increase the signal to noise ratio of the data. The linewidths of the 13 CO(2-1) spectra were corrected for the broadening due to optical depth effects (see explanations in Sect. 4). The median value of σ tot is 0.25 km s -1 with a narrow dispersion of 0.03 km s -1 .

Fig. 8 .

 8 Fig. 8. Effective Jeans diameter (D J,eff = 2σ 2tot /GΣ 0 ) versus filament intrinsic central column density. The vertical dashed line and the grey band mark the border zone between subcritical and supercritical filaments (cf. Fig.6). The dotted line running from the top left to the bottom right corresponds to the thermal Jeans diameter (2R J ). The blue solid lines correspond to the best fits for the subcritical (D sub J,eff ∝ N H 2 -1.1 ± 0.1 )

Fig. A. 1 .

 1 Fig. A.1. Total velocity dispersion versus observed central column density. Similar plot to Fig. 6 for the filaments observed in IC 5146, Aquila and Polaris. The six C 18 O spectra showing two or three velocity components are flagged in yellow.

Fig. A. 5 .

 5 Fig. A.5. Curvelet component of the column density map of the Aquila field around the Serpens South filaments taken fromAndré et al. (2010). The curvelet component is obtained using a morphological component analysis algorithm (MCA, fromStarck et al. 2003) which enhances the contrast of the filamentary structure against the more diffuse background of the cloud. The positions of the observed spectra are plotted in red squares and blue triangles for N 2 H + and C 18 O, respectively. The numbers correspond to the filaments listed in Table1and are given at the upper right of the spectra in Fig. A.4 and some of the spectra of Fig. A.7.

  Fig. A.5. Curvelet component of the column density map of the Aquila field around the Serpens South filaments taken fromAndré et al. (2010). The curvelet component is obtained using a morphological component analysis algorithm (MCA, fromStarck et al. 2003) which enhances the contrast of the filamentary structure against the more diffuse background of the cloud. The positions of the observed spectra are plotted in red squares and blue triangles for N 2 H + and C 18 O, respectively. The numbers correspond to the filaments listed in Table1and are given at the upper right of the spectra in Fig. A.4 and some of the spectra of Fig. A.7.

Fig. A. 6 .

 6 Fig. A.6. Continuation of Fig. A.4 presenting the N 2 H + (1-0) spectra observed toward 9 filaments in Aquila.

Fig

  Fig. A.9. N 2 H + (1-0) spectrum observed toward filament 1 in Polaris. The corresponding hyperfine structure Gaussian fit is highlighted in blue.

Fig

  Fig. A.10. Curvelet component of the column density map of a subregion in the Polaris Cloud taken from André et al. (2010). The positions of the observed spectra are plotted in red squares and blue triangles for N 2 H + and C 18 O, respectively. The numbers correspond to the filaments listed in Table1.

Fig. A. 11 .

 11 Fig. A.11. C 18 O(1-0) spectra observed toward 6 filaments in Polaris.The corresponding Gaussian fits are highlighted in blue. Note 1: the spectrum observed toward filament 4 shows two velocity components. Both velocity components have similar linewidths (0.19 km s -1 and 0.14 km s -1 , respectively), they are separated by ∼0.8 km s -1 . They most probably belong to two filaments which are seen on the same line-of-sight (but a large map would be needed to study the kinematics of the field). The component observed at -4.5 km s -1 has been taken to be representative of filament 4.

Table 1 .

 1 Summary of derived parameters of the observed spectra toward the filaments of our sample.

	Field/tracer Filament	Ra	Dec	N H 2	M line	V LSR	T peak	rms	σ obs	σ tot	σ err tot

  Columns 3 and 4: Equatorial coordinates providing the positions of the observed spectra. Columns 5 and 6: Projected column density and mass per unit length for each filament measured from Herschel column density maps. The listed values have not been corrected by the 60% systematic overestimation of true column densities and masses per unit length. All the filaments detected in N 2 H + are supercritical while most of the filaments observed in C 18 O(1-0) and (2-1) and not detected in N 2 H + are subcritical. Columns 7, 8 and 10: parameters of single and multiple Gaussian fits to the C 18 Oa n dN 2 H + observed spectra, respectively. For the N 2 H + spectra T peak (Col. 8) corresponds to the peak temperature of the isolated component of the N 2 H + multiplet. Column 9 gives the corresponding rms. The reference positions used for the sposition-switched N 2 H + observations were selected from the Herschel images and correspond to 21:

								9)	(10)	(11) (12)
		1	18:29:39.5 -01:50:14	213.4	252	7.70	1.43 0.67	0.60	0.64 0.01
		2	18:32:38.5 -02:50:24	40.9	176	7.85	0.67 0.22	0.30	0.36 0.01
		3	18:33:15.2 -02:48:14	13.0	35	6.42	0.65 0.13	0.11	0.23 0.01
		4	18:33:06.7 -02:47:50	11.6	32	7.04	0.26 0.08	0.11	0.23 0.01
		5	18:33:22.5 -02:42:30	12.8	20	6.84	0.16 0.05	0.15	0.25 0.01
		6	18:29:36.6 -01:51:55	13.2	39	7.63	0.12 0.05	0.11	0.23 0.01
		7	18:29:50.4 -01:54:47	37.9	123	7.67	1.18 0.25	0.18	0.27 0.01
	Aquila	8	18:30:27.6 -01:55:32	35.2	139	7.42	0.41 0.13	0.20	0.28 0.01
	N 2 H + (1-0)	9	18:30:19.1 -01:53:09	33.9	90	7.72	0.18 0.11	0.41	0.45 0.01
		10	18:29:08.9 -01:42:40	29.9	121	7.05	0.76 0.18	0.21	0.29 0.01
		11	18:29:19.9 -01:52:25	25.1	68	7.27	0.32 0.09	0.23	0.31 0.01
		12	18:30:03.9 -02:04:29	20.3	73	7.12	0.68 0.19	0.19	0.27 0.06
		13	18:29:55.8 -02:00:18	74.5	182	7.09	1.39 0.60	0.42	0.47 0.01
		14	18:29:58.1 -01:57:38	55.8	265	7.41	0.45 0.21	0.44	0.48 0.01
		15	18:28:50.1 -01:34:52	114.1	288	7.90	0.79 0.23	0.36	0.41 0.01
		16	18:30:03.5 -02:02:26	10.9	35	7.36	0.37 0.07	0.10	0.22 0.01
		17	18:29:55.4 -01:37:02	12.0	40	8.29	0.53 0.12	0.14	0.24 0.01
		18	18:29:34.6 -01:55:15	9.2	38	7.34	0.10 0.04	0.15	0.25 0.01
		2	21:46:28.1	47:17:34	25.3	80	4.69	0.61 0.17	0.19	0.27 0.01
	IC 5146	6	21:47:17.2	47:33:01	17.0	109	3.69	0.32 0.11	0.25	0.32 0.01
	N 2 H + (1-0)	12	21:44:58.4	47:34:15	14.7	41	4.11	0.40 0.13	0.19	0.28 0.01
		13	21:47:07.9	47:39:30	17.1	50	4.35	0.25 0.10	0.22	0.30 0.01
	Polaris	1	02:00:29.0	87:41:54	8.0	16	-4.62	0.18 0.01	0.13	0.24 0.01
	N 2 H + (1-0)
		19	18:29:53.5 -03:43:39	7.9	10	4.88	1.22 0.06	0.23	0.31 0.01
		20	18:28:58.2 -01:22:51	7.8	28	7.17	1.10 0.04	0.30	0.36 0.01
		21	18:31:06.4 -01:18:58	9.9	20	7.77	1.10 0.05	0.26	0.33 0.01
	Aquila	22	18:33:15.8 -02:47:15	7.2	51	6.47	1.26 0.03	0.13	0.24 0.01
	C 18 O(1-0)	23	18:22:51.1 -03:26:38	6.3	12	5.78	0.36 0.03	0.28	0.34 0.01
		24	18:23:31.3 -03:03:43	7.2	15	6.87	0.89 0.03	0.27	0.34 0.01
		25	18:29:45.6 -02:32:21	7.3	22	6.52	0.94 0.03	0.15	0.25 0.01
		26	18:32:53.1 -01:29:42	7.0	21	7.28	0.28 0.03	0.16	0.26 0.01
		3	21:46:42.1	47:44:26	2.8	3	1.89	0.80 0.08	0.17	0.26 0.01
		5	21:45:15.3	47:46:52	3.0	11	1.42	0.63 0.09	0.24	0.31 0.02
	IC 5146	7	21:44:37.1	47:42:39	5.2	11	3.71	0.76 0.09	0.29	0.35 0.02
	C 18 O(2-1)	10	21:43:29.1	47:19:29	4.3	16	4.21	0.79 0.06	0.13	0.24 0.01
		11	21:43:30.9	47:47:12	3.1	9	3.29	0.42 0.08	0.11	0.23 0.01
		14	21:47:48.1	48:04:53	3.9	15	0.24	0.89 0.07	0.12	0.23 0.01
		25	21:44:18.8	47:32:25	5.7	17	2.32	0.47 0.09	0.27	0.34 0.03
		2	01:58:59.7	87:39:33	9.6	16	-4.29	2.65 0.04	0.22	0.29 0.01
		3	01:40:31.4	87:46:19	2.4	15	-4.14	0.28 0.04	0.22	0.30 0.02
	Polaris	4	01:31:21.7	87:45:56	2.7	25	-4.48	0.27 0.03	0.19	0.27 0.02
	C 18 O(1-0)	5	02:00:49.1	87:35:18	3.4	15	-4.08	0.62 0.04	0.17	0.30 0.01
		6	01:50:48.0	87:43:47	1.5	3	-3.46	0.19 0.04	0.19	0.28 0.02
		7	03:36:47.4	88:08:13	2.3	16	-3.39	0.12 0.05	0.08	0.26 0.01
	Notes. 47:22.6 47:45:59 for IC 5146, 01:56:30.26 87:43:28 for Polaris,
	18:28:23.9 -02:00:22 for Aquila, except for filaments 3, 4 and 5 for which a reference position at 18:33:27.75 -02:51:22 was used. Columns 10 and 11: σ obs and σ tot are derived as explained in Sect. 3.

http://www.iram.fr/IRAMFR/GILDAS

Due to the presence of apparent intersections of filamentary structures in the plane of the sky, there is sometimes some ambiguity in the definition of an "individual filament".

The column densities were derived consistently for all filaments observed with Herschel in Polaris, IC 5146, Aquila, and DR21. The same dust opacity law was used to derive the column density of NGC 2264C from MAMBO 1.2 mm continuum observations. The column density estimates for NGC 2264C and DR21 are more uncertain due to the presence of prominent cores along the filament crests (cf.Peretto et al. 2006;Motte et al. 2007).

Pages 11 to 15 are available in the electronic edition of the journal at http://www.aanda.org

Acknowledgements. D.A. wishes to thank the staff at the IRAM 30 m telescope at Pico Veleta for the support during the observations and Q. Nguy ẽn Lu'o'ng for his help in reducing the IRAM 30 m data using GILDAS/CLASS. This work has benefited from the support of the European Research Council under the European Union's Seventh Framework Programme (Grant Agreement No. 291294) and of the French National Research Agency (Grant No. ANR-11-BS56-0010).

Appendix A: Velocity components of the observed filaments

The N 2 H + (1-0) and C 18 O(1-0)/(2-1) spectra observed toward the sample of filaments studied in this work are shown in Appendix A. The 23 filaments detected in N 2 H + (1-0) in Aquila, IC 5146, and Polaris all show a single velocity component (cf. Figs. A.2, A.4, A.6 and A.9). This is consistent with the view that the filaments are velocity coherent structures (e.g., Hacar & Tafalla 2011). While most of the filaments detected in C 18 O show one velocity component, some of them show multiple (two or three) C 18 O velocity components (such as filaments 20,22,23,25,26 in Aquila and filament 4 in Polarissee Figs. A.7,A.8 and A.11). Some of these multiple C 18 Ov elocity components may be due to line-of-sight mixing of emissions not necessarily tracing the selected filament, but background and/or foreground structures (especially in the case of the Aquila Rift which is a region close to the Galactic plane).

Assuming that interstellar filaments are velocity-coherent structures, the most likely velocity component associated with each filament was considered (selecting the velocity component with the velocity closest to that of neighboring filaments). Its properties were then used in the discussion of the results presented in this paper (see detailed explanations in the captions of the figures presenting the observed spectra).

In order to investigate how the uncertainty in selecting the relevant velocity component affects our results, the 6 filaments with spectra showing more than one C 18 O velocity component were flagged in the sample. 1. The corresponding single-component hyperfine structure Gaussian fits are highlighted in blue. Note that the y-axis scale is not the same for all spectra, but has been adjusted to match the peak temperatures which vary between ∼0.2 to >3K. The component at V LSR = 7.17 km s -1 have been chosen to be tracing the filament, since this component has a velocity component which is closer to the average velocity (V LSR = 7.4k ms -1 ) of the filaments in Aquila which show one velocity component. Note 2: the spectrum observed at the position of filament 22 shows two velocity components, one centered at 6.5 km s -1 and the other at 7.7 km s -1 . The strongest component has been selected for similar reasons to that explained for filament 20. Note 3: the spectrum observed at the position of filament 25 has two velocity components (at 5.0 km s -1 and 6.52 km s -1 , respectively). The strongest component has been selected for similar reasons to that explained for filament 20. Note 4: the spectrum observed at the position of filament 26 shows three velocity components centered at 7.28 km s -1 ,8 .4 5k ms -1 and9.17kms -1 , respectively. The component centered at 7.28 km s -1 has been selected because it is the component with the closest V LSR value to the mean V LSR of neighboring filaments in Aquila.


Thomas Leroy
email: thomas.leroy@ljll.math.upmc.fr

A well-balanced scheme for a transport equation with varying velocity and irregular opacities

An equation for photons containing a frequency drift term modeling the Doppler effects and an emission absorption coefficient is considered. On the one hand, this frequency drift term involves a coefficient κ relative to the fluid velocity with which the photons interact. This coefficient may vanish, leading to numerical difficulties. On the second hand, the emission absorption σ is known to be very irregular with respect to the frequency, and thus the design of well-balanced scheme is a mathematical issue. In this work are presented new numerical results for the spectrally well-balanced scheme in the context of highly irregular opacities.

Introduction

The aim of this work is to extend the recent study [START_REF] Leroy | Després A well-balanced scheme for a transport equation with varying velocity arising in relativistic transfer equation[END_REF] to highly irregular emission absorption coefficients. The coupling between the transport equation for photons and an equation describing the evolution of a fluid moving at velocity u involves several parameters. In particular, the emission absorption coefficient σ is known to be very irregular with respect to the frequency. It yields mathematical difficulties, in particular to prove theoretical results as existence of solutions (see for example [START_REF] Bardos | The Rosseland Approximation for the Radiative Transfer Equations Communications on Pure and Applied Mathematics[END_REF][START_REF] Bardos | The Nonaccretive Radiative Transfer Equations: Existence of Solutions and Rosseland Approximation Journal of Functional Analysis[END_REF]). This kind of equations is complicated to solve, and thus reduced models is an issue. Recently ( [START_REF] Leroy | Relativistic transfer equations: maximum principle and convergence to the nonequilibrium regime Work in progress[END_REF]), the nonequilibrium regime in the context of a relativistic gas of photons interacting with a fluid of velocity u has been studied, and the convergence in this regime of the solution of the relativistic transfer equation to the solution of a diffusion equation involving a frequency drift term has been proved. The corresponding homogeneous (in space) equation is

   ∂ t ρ = κ 3 ν∂ ν ρ + σ(ν)(B(ν) -ρ) in [0, T ] × R + ν , ρ(0, ν) = ρ in (ν), (1.1)
where ρ = ρ(t, ν) is the density of photons at time t and at frequency ν, B(ν) = ν 3 e ν -1 -1 is the Planck function, σ is the emission absorption coefficient and κ ∈ [-κ * , κ * ], κ * > 0 is the divergence of the fluid velocity. The design of costless numerical scheme for this equation is an interesting issue. Since radiative phenomena are very fast compared to hydrodynamic ones, wellbalanced schemes are a possible solution. As mention in [START_REF] Leroy | Després A well-balanced scheme for a transport equation with varying velocity arising in relativistic transfer equation[END_REF], classical "Greenberg-Leroux" type well-balanced schemes [START_REF] Greenberg | A Well-Balanced Scheme for the Numerical Processing of Source Terms in Hyperbolic Equations[END_REF] (see also [START_REF] Gosse | An asymptotic preserving well-balanced scheme for the hyperbolic heat equations C[END_REF] for a recent review on the topic) constructed by discretizing the source terms at the interfaces and by using a Godunov scheme on the obtained equation are not consistent with equation (1.1) in the regime κ → 0. This explains the introduction of the spectrally well-balanced scheme (SWB), for which a uniform (with respect to the waves velocity κ) convergence result has been provided. On the other hand and as mentioned before, the emission absorption coefficient is highly irregular, and plays an important role from the physical point of view. It is thus necessary to report faithfully this behavior. In a first approach a very fine mesh could be used, but it is not possible in general in photons' transport. Indeed the distribution function of the photons is often rather peaked in frequency. In this context and due to the frequency shift, the classical upwind scheme is not accurate enough.

In this paper the behavior of the spectrally well-balanced scheme is studied for several irregular opacities, and its numerical solution is compared with the Greenberg Leroux and upwind schemes. In particular numerical tests are performed with a coefficient σ constructed as the sum of a regular part and dirac masses, which seems to be realistic, and with the so called Cramer's opacity. The well-balanced property and the consistency as κ → 0 ± of the spectrally well-balanced scheme are shown to be preserved for such opacities. This paper is organized as follow. In a first part the different schemes of interest and their properties are recalled. In particular the behavior of each schemes as κ → 0 and their dependence (via the CFL) with respect to the emission absorption coefficient σ are pointed out. In a second part numerical results for different configurations of the parameters κ and σ, involving in particular irregular opacities, are presented. The consistency of the SWB scheme as κ → 0 even for irregular σ is shown numerically. Finally, the interest of the SWB scheme in relation to the upwind scheme will be illustrated with test cases in which the relative L 1 error of the upwind scheme can reach 35%.

A review of well-balanced schemes

Let us consider a frequency domain D = [ε, ν], for given 0 < ε < ν < +∞ and, for 1 ≤ j ≤ N , an irregular mesh defined by (N + 1) points 0 = ν 1 2 < ... < ν N + 1 2 = ν * . The middle of the j-th frequency band is denoted as ν j , i.e. ν j = (ν j-1 2 + ν j+ 1 2 )/2 and its length is denoted as ∆ν j . We also define the dual (j + 1 2 )-th frequency band as the cell [ν j , ν j+1 ], which length is denoted ∆ν j+ 1 2 and h = max j ∆ν j . The existence of a constant C such that ∀j ∈ {1, ..N } , 0 < Ch ≤ ∆ν j is assumed. Since the aim of this paper is the design of well-balanced schemes, it is necessary to study the stationary equation associated to (1.1)

   κ 3 ν∂ ν ρ + σ(ν)(B(ν) -ρ) = 0 in [0, ν], ρ(ν * ) = ρ * , (2.1)
where ν * = ν if κ > 0, which yields a transport of the photons toward the frequency 0, and ν * = 0 if κ < 0, which yields a transport of the photons toward the frequency ν. Obviously this equation can be solved and one has

ρ ν; ρ * , ν * = ρ * e -3 κ ν * ν σ(s) s ds + 3 κ ν * ν σ(s)B(s) s e -3 κ s ν σ(τ ) τ dτ ds. (2.2)
This solution as a limit as κ → 0 ± :

lim κ→0 ± ρ ν; ρ * , ν * = B(ν). (2.3)
This stationary solution will be used in the next parts to construct well-balanced schemes.

" Classical " well-balanced scheme

In this part a first well-balanced scheme for equation (1.1) is introduced. Its construction uses the method developed by J.M. Greenberg and A.Y. Leroux [START_REF] Greenberg | A Well-Balanced Scheme for the Numerical Processing of Source Terms in Hyperbolic Equations[END_REF]: discretizing the source terms at the interfaces, one uses a Godunov scheme on the obtained equation. It yields the following scheme, denoted as WB1 scheme

       d dt ρ j = κ 3 ν j ρ(ν j ; ρ j+1 ; ν j+1 ) -ρ j ∆ν j , 1 ≤ j ≤ N, κ > 0, d dt ρ j = κ 3 ν j ρ j -ρ(ν j ; ρ j-1 ; ν j-1 ) ∆ν j , 1 ≤ j ≤ N, κ < 0.
(2.4)

Considering a classical Euler discretization of the time derivative, the CFL is given by ∆t ≤ max j 3∆ν j |κ|ν j .

(2.5)

It is interesting to note that the CFL of this scheme does not depend on σ, but depends on the wave velocity κ. A simple Taylor expansion of the exponential functions in the stationary solution shows the consistence of this scheme for any fixed and non zero κ. On the other hand, this scheme is not consistent as κ → 0 ± . Actually, the limit of the stationary solution in this regime (2.3) yields lim

κ→0 ± d dt ρ j = 0, (2.6)
which obviously is not a consistent discretization of equation (1.1) as κ = 0. This explain the introduction of a new class of well-balanced schemes

Spectrally Well-Balanced scheme

In this part the Spectrally Well-Balanced scheme [START_REF] Leroy | Després A well-balanced scheme for a transport equation with varying velocity arising in relativistic transfer equation[END_REF], whose construction use the integrating factor M (ν

* ; ν) = e -3 κ ν * ν σ(s)s -1 ds , is recalled. It is defined as        d dt ρ j = σ(ν j ) 1 -M (ν j+1 ; ν j ) ρ ν j ; ρ j+1 , ν j+1 -ρ j , 1 ≤ j ≤ N, κ > 0, d dt ρ j = σ(ν j ) 1 -M (ν j-1 ; ν j ) ρ ν j ; ρ j-1 , ν j-1 -ρ j , 1 ≤ j ≤ N, κ < 0.
(2.7)

Considering once again a classical Euler discretization of the time derivative, the CFL is given by

       ∆t ≤ max j 1 -M (ν j+1 ; ν j ) σ(ν j ) , κ > 0, ∆t ≤ max j 1 -M (ν j-1 ; ν j ) σ(ν j ) , κ < 0.
(2.8)

On the contrary to the previous scheme, the CFL of the SWB scheme strongly depends on the L ∞ norm of the coefficient σ. Actually, one has lim σ→∞ ∆t = 0.

On the other hand, it is interesting to see that it is the integrating factor that allows to recover a consistent limit scheme as κ → 0 ± . Actually, the limit (2.3) as κ → 0 ± of the stationary solution leads to lim

κ→0 ± d dt ρ j = σ(ν j ) B(ν j ) -ρ j .
(2.9)

One of the main properties of this scheme is the following uniform (with respect to κ) result.

Lemma 2.1 (Uniform convergence). We define ρ h = (ρ j ) 1≤j≤N and the discrete norm

. L 2 d such that ρ h 2 L 2 d = j ∆ν j+ 1 2 ρ 2 j .
Assume that the initial data satisfies ρ in ∈ H 2 (D). Assume moreover that the emission absorption coefficient satisfies σ ∈ W 2,∞ (D) and that there exists a constant σ * > 0 such that ∀ν ∈ D, σ(ν) ≥ σ * . Then the numerical solution of the scheme (2.7) satisfies the following estimate, where the constant

C(T ) is uniform in κ ∈ (0, κ * ]: ρ ex (t) -ρ h (t) L 2 d ≤ C(T )h , 0 < t < T .

Upwind scheme

In all the forthcoming numerical results, the numerical solutions of the well-balanced schemes (2.4) and (2.7) are compared to the following classical upwind scheme

       d dt ρ j = κ 3 ρ j+1 -ρ j ∆ν j ν j+ 1 2 + σ(ν j ) B(ν j ) -ρ j , 1 ≤ j ≤ N, κ > 0, d dt ρ j = κ 3 ρ j -ρ j-1 ∆ν j ν j-1 2 + σ(ν j ) B(ν j ) -ρ j , 1 ≤ j ≤ N, κ < 0.
(2.10)

Considering once again a classical Euler discretization of the time derivative, the CFL is given by

         ∆t ≤ 3 max j ∆ν j κν j+ 1 2 + 3∆ν j σ(ν j ) , κ > 0, ∆t ≤ 3 max j ∆ν j -κν j-1 2 + 3∆ν j σ(ν j ) , κ < 0.
(2.11)

As for the previous SWB scheme, the upwind scheme strongly depends on the L ∞ norm of σ.

Obviously the upwind scheme is consistent in both regimes ∆ν j → 0 and κ → 0 ± .

Numerical results

In this section numerical results are presented for the " classical " well-balanced (2.4) scheme, the Spectrally well-balanced (2.7) scheme and the upwind (2.10) scheme in several configurations of the coefficients κ and σ. This section is divided in two parts. In a first part the long time behavior and the consistency as κ → 0 ± of the schemes of interest are studied for highly irregular emission absorption coefficients σ. The main point is that the SWB scheme exactly capture the numerical solution as time goes on and is consistent as κ → 0 ± . In a second part we study several configurations of the parameters κ and σ in which the relative L 1 error of the upwind scheme is maximized.

Since this study involves a very simple model, the analytical solution of equation (1.1) can be computed, and its expression is given by

ρ(t, ν) = ρ in (νe κ 3 t )e -3 κ νe κ 3 t ν σ(τ )τ -1 dτ + t 0 σ νe κ 3 s B νe κ 3 s e -3 κ νe κ 3 s ν σ(τ )τ -1 dτ ds, (3.1)
Since one of the regime of interest is κ << 1, it is numerically necessary to use a non singular form of this analytical solution. To this end one remarks that the expression under the integral can be integrated by parts. One finds

ρ(t, ν) = B(ν) + ρ in (νe κ 3 t ) -B νe κ 3 t e -3 κ νe κ 3 t ν σ(τ )τ -1 dτ + κ 3 ν t 0 e κ 3 s B ′ νe κ 3 s e -3 κ νe κ 3 s ν σ(τ )τ -1 dτ ds, (3.2)
which is no longer singular in κ. The same procedure is applied to treat the singularity in the stationary solution (2.2). The following relative L 1 error is introduced

ρ j (t) -ρ(ν j , t) ρ(ν j , T f ) . (3.3)
For all the numerical tests a random mesh composed of 50 cells is considered. The domain is D = [ε, 30], where ε varies between 10 -5 and 10 -1 according to the different parts. The initial conditions are taken as ρ in (ν) = B(ν -15) if κ > 0 and ρ in (ν) = 0 if κ < 0 and the boundary conditions are given by the analytical solution (3.2).

Highly irregular σ

Many physical contexts implies highly irregular emission absorption coefficient σ, whose frequency dependence is sometimes regular and sometimes very peaked. The numerical study of such opacities is thus interesting. Numerically, the coefficient σ is constructed as the sum of a regular part and dirac masses (cf figure 1)

σ(ν) = 150 1 + ν + 1≤j≤20 δ j , (3.4)
where the dirac functions are distributed regularly on the mesh. In this part ε is taken equal to 10 -5 , so the domain is [10 -5 , 30].

Figure 1: Irregular emission absorption coefficient

The choice of σ (3.4) yields an exact numerical integration of all the integrals in the analytical solution (3.2). In the same way, the integrals in the expression of the stationary solutions (2.2) are numerically computed exactly. The first property to verify numerically is the well-balanced property of the SWB (and WB1) scheme. In figure 2 are displayed the time evolution of the L 1 relative error of the SWB, WB1 and upwind scheme for κ = 1 and κ = -1. As expected, both the SWB and the WB1 schemes exactly capture the stationary solution, but it is interesting to notice that the WB1 scheme is much slower than the SWB scheme. As explain in the previous section, the SWB scheme is uniformly (in κ) convergent, and this property has to be checked numerically for this choice of σ. Figure 3 displays the evolution, as κ → 0, of the L 1 error between the solutions of the SWB, WB1 and the upwind schemes and the numerical solution of the following scheme, consistent with equation (1.1) as κ = 0:

d dt ρ j = σ(ν j ) B(ν j ) -ρ j .
(3.5)

Cramer's opacity

In this part we study the behavior of the schemes with the Cramer's opacity:

σ(ν) = 1 -e -ν ν 3 . (3.6)
This opacity leads to σ(ν)B(ν) = e -ν , and has a singularity as ν → 0. For this coefficient the analytical value of the integrals in the expression of solutions (1.1) and (2.1) can no longer be computed. Numerically, they are computed by classical five points Gaussian quadratures. In this part ε is taken equal to 10 -1 and thus the domain is [10 -1 ; 30]. The same study than for the previous opacity is performed, and as before the long time behavior of the numerical solutions of the schemes is studied. In figure 5 are displayed the time evolution of the relative L 1 error of the schemes. Once again and as expected, the SWB and WB1 schemes exactly capture the stationary solution, while the relative error of the upwind scheme is 18% for κ = 1 and 7% for κ = -1. As before, the behavior of the schemes for κ → 0 ± is studied. Figure 6 displays the error between the schemes and the numerical solution of the scheme (3.5) as κ → 0 ± . It shows that the SWB scheme is consistent with equation (1.1), even for small coefficient κ.

SWB vs upwind

In this part the advantage of the SWB scheme compared with the upwind scheme is highlighted. Indeed, for certain value of κ and σ, the upwind L 1 error can become significant. In the previous part, for σ defined in (3.4) and κ = ±1, the upwind relative L 1 error was of order 1%. On the other hand we have seen that for the Cramer's opacity and for κ = 1, the relative L 1 error of the upwind scheme reaches 18%, which is non negligible. Two test cases are presented here in which the upwind relative L 1 error is significant. In both the case we take σ as a piecewise constant function. In the first following one, the discontinuity is placed at the middle of the mesh (figure 7). For the second test case the emission absorption coefficient is still taken as a piecewise constant function but we took a negative κ. The initial condition are such that ρ j (0) = 0. In figure 10 are displayed the evolution of the relative L 1 error as time goes on, for several value of κ. As previously, we see that the relative L 1 error of the upwind scheme increases with κ, and reaches roughly 20% for κ = -10.

Conclusion

The long time behavior of the Spectrally Well-Balanced scheme has be studied for several irregular opacities, including highly peaked emission absorption coefficient and the Cramer's opacity. In particular the well-balanced property of the scheme was shown to be preserved, and its consistency as the velocity speed κ tends to 0 was also studied numerically for these opacities, and the numerical results agree with the theoretical result (theorem 2.1). The case κ < 0 has also been studied. Finally, the advantage of the SWB scheme with respect to the upwind scheme has been shown in two test cases for different values of the coefficients κ and σ.

Figure 2 :

 2 Figure 2: L 1 Relative error vs time. Top: κ = 1, bottom: κ = -1

Figure 3 :

 3 Figure 3: L 1 error vs κ. Left: κ > 0, t=0.1. Right: κ < 0, t=0.02

Figure 4 :

 4 Figure 4: L 1 error vs N , (Log,Log) scale. Left: κ > 0, t=0.1. Right: κ < 0, t=0.02

Figure 5 :

 5 Figure 5: Relative L 1 error vs time. Left: κ = 1, right: κ = -1

Figure 6 :

 6 Figure 6: L 1 error vs κ. Left: κ > 0, t=0.5. Right: κ < 0, t=2

Figure 7 :

 7 Figure 7: Piecewise constant emission absorption coefficient

Figure 8 :

 8 Figure 8: Relative L 1 error vs time

Figure 9 :

 9 Figure 9: Corresponding solutions at t=30

Figure 10 :

 10 Figure 10: Relative L 1 error vs time

Figure 11 :

 11 Figure 11: Corresponding solutions at t=30