© ESO, 2015

1. Introduction

It is widely accepted that the gas-to-dust mass ratio (g2d) in our Galaxy has a value of about 100. Of course, this is an average value and it does, in particular, pertain to the diffuse interstellar medium (ISM). However, this value is also commonly adopted for dense cloud environments, such as star-forming regions, and the methods used to determine that parameter are described for instance by Kenyon et al. (1998), Draine (2003), and Vuong et al. (2003). Based on observations of CO isotopologues and stellar extinction values in the nearby star-forming regions in Taurus and Ophiuchus, Frerking et al. (1982) established the still widely used molecular gas-dust relations, that is, N(^XC^YO)∝A_V, A_V ≲ 20 mag, and ¹²C¹⁶O/H₂ ~ 8.5 × 10^-5. For higher A_V, saturation makes these relations less useful.

Young & Scoville (1991) discussed the high average g2d of ~600 obtained for molecular clouds, where gas masses refer to H I + H₂. The latter was derived from CO observations and the dust masses were determined from IRAS observations. Furthermore, based on data for H I and H II regions, Spitzer (1978) quoted the rather wide range of 20 to 700 for this parameter. A value of about 450 was found by Liseau et al. (1995) for L 1688 using CS maps and 60 μm IRAS data. A high g2d would also be consistent with the low oxygen abundance in ρ Oph A found by Liseau & Justtanont (2009). Figure 11 of Brinchmann et al. (2013), which is based on a large number of galaxies, supports this view on larger, extragalactic scales.

One of the reasons for these high g2d values could be that large portions of the dust are at temperatures considerably lower than 30 K and, hence, remained undetected by IRAS (Young & Scoville 1991). However, at temperatures below about 25 K, the abundance of the dominant gas tracer, CO and its isotopologues, could also become highly diminished by freeze-out (e.g., Lippok et al. 2013, and references therein). One might therefore select spectral lines of a molecule that is seen in regions where CO is thought to be depleted. One such molecule is diazenylium, N₂H⁺ (e.g., Caselli et al. 1999), which is formed from N₂ (which remains in the gas phase to somewhat lower temperatures than CO) and is removed from the gas phase mainly by CO (Bergin et al. 2002; Öberg et al. 2005). However, as pointed out by Bergin et al. (2002), even N₂H⁺ (and N₂) may become depleted and thus leave very few molecular probes (perhaps only H₂D⁺, Brünken et al. 2014) to study the centres of the densest cores.

The rotational transitions (J = 3−2) and (J = 6−5) of N₂H⁺, at 280 GHz and 560 GHz and with upper level energies of 27 K and 94 K, respectively, could fulfil these requirements. In particular, the (1−0) line at 93 GHz, having an upper level energy similar to that of CO, has been widely used for dense interstellar clouds, including the ρ Oph A region (e.g., Di Francesco et al. 2004). However, E₁/k is only 4.5 K and at considerably higher temperatures, the population of the J = 0 and J = 1 levels changes only slowly and does as such not contribute much to the understanding of the overall excitation of the molecule.

Friesen et al. (2014) have used the Atacama Large Millimeter/submillimeter Array (ALMA) to observe ρ Oph A in the N₂H⁺ (4−3) lines toward a small region (< 20′′ × 20′′), but the paper focused on the continuum data. There are no line profiles nor line maps shown for these data, which, according to Friesen et al. (2014), were either self-absorbed or resolved out by the interferometer. The authors also noted that line and continuum sources did not coincide, neither for SM 1 nor for SM 1N, a point of potential importance (see below, Sect. 4.4). The primary beam of these observations is smaller than the single-dish beams of our observations.

As revealed by pure-rotational H₂ line emission, relatively high temperatures, ≳10³ K, are present in the interface regions toward the dominating radiation source S 1 (Fig. 1). In addition, Liseau et al. (1999) found extended [O i] 63 μm, [O i] 145 μm and [C ii] 157 μm emission, which is indicative of temperatures in excess of 100 K. These phenomena probably refer to the outer layers of the cloud. The N₂H⁺ (6–5) line could be used to probe that temperature regime in a coherent map.

Already early on, the actively star-forming ρ Oph clouds attracted wide attention, largely because the nearest stellar cluster is forming there. This is the dominant star formation mode, in contrast to more isolated star formation as in the Taurus-Auriga clouds, for example, that is contributing only little to the stellar initial mass function (IMF) of the Galaxy. Not surprisingly therefore, a huge body of literature exists, and we refer to the review by Wilking et al. (2008) and also to a few recent specific papers on the distribution of molecular gas and icy grains (Shuping et al. 2000; Pontoppidan 2006; André et al. 2007; White et al. 2015). For these clouds in ρ Oph, the visual extinction can locally reach extremely high values (maxA_V> 10² mag).

The nearby 3′ (0.1 pc) ρ Oph A core (Loren et al. 1990) is unique in the sense that all phases of the earliest phases of stellar evolution are present in a very compact region (see Figs. 1 and 2) – from starless clumps, Class −1 (SM 1), via Class 0 (VLA 1623), Class I (GSS 30-1) and Class II (LFAM 3) to Class  III (GSS 32) (see, e.g., Comeron et al. 1993; Bontemps et al. 2001; White et al. 2015). Potentially, this fact may make it possible to follow the time evolution of the dust in various locations in the core, including grain growth and changes in composition, such as ice coatings onto initially “bare” grains. Such a scenario could have considerable effect on theories of planetary formation (e.g., the location of the snow line in the circumstellar disks) and help to relax tight timescale requirements, particularly in the outer, less dense regions of the disks. The status of the current understanding of the dust evolution in disks has recently been summarized by Testi et al. (2014).

In this paper, we focus on the gas and dust contents of the nearby star-forming ρ Oph A core, observed in the far-infrared (FIR) and sub-millimeter spectral regions (submm). For the distance we adopt the accurately determined value with the Very Long Baseline Array (VLBA) by Loinard et al. (2008, viz. . The organization is as follows: in Sect. 2, our mapping observations of freeze-out insensitive molecules are described. These observations have been obtained from space with Herschel (Pilbratt et al. 2010) and from the ground with APEX at 5100 m altitude. There, maps in the submm continuum have also been obtained. This section also describes the reduction of these data, and the results are presented in Sect. 3. These are discussed in Sect. 4, illustrating specifically the aspects of spectral line overlap and also of geometry in the radiative transfer problems. We compare the results obtained with tools widely in use with those arrived at with our own developments. The conclusions follow in Sect. 5.

2. Observations and data reduction

2.1. SHeFI-APEX-2 and Herschel-HIFI

2.1.1. APEX-2

The Atacama Pathfinder Experiment (APEX) is a 12 m single-dish telescope situated at 5100 m altitude in Chile. The pointing accuracy of the telescope has been determined at 2′′ (rms). We used the Swedish heterodyne facility instrument APEX-2 (Vassilev et al. 2008) for the raster-mapping observations of ρ Oph A in the J = 3−2 rotational transition of the N₂H⁺ molecule. The central rest frequency is at 279 515.06711 MHz and the average zenith optical depth was 0.04. At 279.5 GHz, the beam size is FWHM = 22′′, the main beam efficiency is η_mb = 0.73, and the LSB receiver temperature was 125 K. We used the Fast Fourier Transform Spectrometer with 73.6 kHz wide channels, resulting in a velocity resolution of 0.082 km s^-1 per channel and a spectral coverage corresponding to 200 km s^-1. The data reduction involved the fitting and removal of the baselines and was made with the locally available software package xs¹.

2.1.2. HIFI

The (J = 6−5) line of N₂H⁺ was observed at the same time as our Herschel H₂O maps (2012, GT2_rliseau_1). The observational details regarding the partial map data of 2010 (Table 1) that were obtained within the KPGT_evandish_1 programme were already described in detail by Bjerkeli et al. (2012). In their Fig. 1, offsets are relative to the (0, 0)-position at α₂₀₀₀ = 16^h26^m264, δ₂₀₀₀ = −24°24′31′′. The more extended mapping observations in 2012 were made on-the-fly, with a reference position 14′ away. The data reduction was made within HIPE 10². For the final composite map, the zero-offset is at α₂₀₀₀ = 16^h26^m2778, δ₂₀₀₀ = −24°23′503, which is slightly offset from the APEX origin by (−16, −67) and which is smaller than the beam size at 559 GHz, that is, 38′′. At this frequency in HIFI-band 1, the main beam efficiency³ is η_mb = 0.63.

2.2. LABOCA and SABOCA

The continuum observations were carried out with two bolometer multi-pixel cameras, viz. the Large Apex Bolometer Camera LABOCA (Siringo et al. 2009) and the Submillimetre Array, that is, SABOCA (Siringo et al. 2010). The effective operating frequencies correspond to 870 μm and 350 μm, with about 150 and 50 μm wide filters, respectively⁴. The resolutions are 195 and 75, respectively. At the two wavelengths, the maps of ρ Oph A span almost 40′× 40′ and 8′× 8′, respectively (Fig. 3), and in both maps the centre coordinates are at RA = 16^h26^m2790 and Dec = −24°23′570 for the epoch of J2000.0.

The 114 circular field of view (FOV) of LABOCA is determined by the 295 bolometers, which are arranged with 36′′ spacings in between them. The number of loosely spaced bolometers of SABOCA is 39 and its FOV is 15. The FOV is under-sampled for both cameras and is filled by repetitive continuous scanning of the mapped region. For both LABOCA and SABOCA, the data reductions were made with CRUSH⁵. CRUSH uses the actual atmospheric opacities at the time of observation. For the two data sets, average zenith optical depths were 0.42 at 870 μm and 1.07 at 350 μm, with a range of 0.55 to 1.20 for the latter.

3. Results

3.1. Spectral line data

The maps of main-beam corrected integrated intensity, , of N₂H⁺ (6−5) and (3−2) are shown in Fig. 2. The mapping of the (6−5) line with HIFI was made at an angle with respect to the equatorial coordinate system, so that the map had to be rotated to be properly aligned with the APEX data (Fig. 2c). In Figs. 2a and b, the shown pixel sizes correspond to the respective samplings, and in b, the superposed contours show the (3−2) data after convolution to the (6−5) resolution (38′′).

3.2. Continuum data

The reduced data are shown in Fig. 3 for both LABOCA and SABOCA, where intensity units are Jy beam^-1. These intensity scales assume Gaussian beams given in Table 1, and in the figure, the contour levels for LABOCA are (1) –1.21, (2) 0.13, (3) 0.59, (4) 1.20, (5) 2.02, (6) 3.13, (7) 4.61, (8) 6.60 Jy beam^-1. For SABOCA, the corresponding contours are (1) –1.01, (2) 0.31, (3) 0.74, (4) 1.31, (5) 2.06, (6) 3.06, (7) 4.38, (8) 6.12, (9) 8.44, (10) 11.50, (11) 15.55, (12) 20.92, (13) 28.03 Jy beam^-1. The value of 1 Jy beam^-1 corresponds to an intensity of about 100 MJy sr^-1 for LABOCA and 600 MJy sr^-1 for SABOCA, and with 1σ of the rms noise corresponding to 5 and 30 MJy sr^-1, respectively.

At the higher resolution of the SABOCA observations, individual sources can be identified, and their equatorial coordinates are presented in Table 2. The rms pointing accuracy of the APEX telescope is about 2′′. The observed field at 350 μm is much smaller than that at 870 μm (Fig. 3). For quantitative applications, the SABOCA data have to be degraded to the appropriate LABOCA values, that is, 195, and re-sampled onto a 4′′ grid. This is illustrated in the middle panels of Fig. 3, where SABOCA contours are shown superposed on the LABOCA intensities. In addition, the IRAM 1.3 mm data of Motte et al. (1998) are shown, also convolved to the 195 resolution of LABOCA and sampled onto the 4′′ regular grid. For the comparison with the APEX results, the map of ρ Oph A was shifted to the same centre coordinates.

4. Discussion

We start our discussion with determining the physical parameters for the gas, in particular T_kin and N(H₂), using spectral line maps of (presumably) optically thin emission, and C¹⁸O has often been used to guarantee this. However, as shown by Liseau et al. (2010) on the basis of C¹⁸O (3−2) and ¹³C¹⁸O (3−2) observations of ρ Oph A, optical depths do exceed unity even in the rare isotopologues of CO. In those cases, the observed radiation does not carry information about the entire core but is restricted to layers closer to the surface, leading to erroneous estimates of the overall temperature and density distribution. For instance, temperatures based on such observations are around the (for a dense core) surprisingly high value of ~25 K. At temperatures below this, which are commonly encountered in dense cores, the CO gas will be depleted because of freeze-out of the molecules onto the surfaces of the accompanying grains. In contrast, N₂H⁺ molecules remain in the gas phase down to much lower temperatures (Lippok et al. 2013, and references therein). Except for a few exceptions, only the (1−0) line and its hyperfine components have been observed and the analysis generally assumed thermodynamic equilibrium (LTE). Here, we make multi-transition analyses and treat the radiative transfer in non-LTE. The implications for our observations are described in the next section.

4.1. Gas temperature, opacity, and mass

For large numerical grids, the publicly available computer code radex (van der Tak et al. 2007) is very suitable. In its present form, the code is unable to properly handle the transfer in optically thick lines with extremely closely spaced transitions, however (the “line overlap” problem, Fig. 6). In addition, by construction, the code will not return the proper shapes of the lines. To surmount these difficulties, we used for the analysis an accelerated lambda iteration (ALI) code that takes line-overlap explicitly into account and also correctly computes the line profiles with which the observations can be compared directly.

The upper level energy of the N₂H⁺ (3−2) transition corresponds to nearly 28 K, close to previous estimates of the kinetic gas temperature. The observed (3−2) lines exhibit the remarkable feature that the satellite lines reveal strengths that are almost similar to that of the main line (Fig. 4), that is, much stronger than expected for optically thin radiation in LTE (Fig. 6). In fact, this strongly suggests optically very thick emission in the overlapping components of the central, spectrally non-resolved, main line. This line anomaly is in contrast to the observed (6–5) line profile, which appears nearly Gaussian (however, the spectral resolution is lower). The J = 6 level is about 94 K above ground.

The Onsala ALI code (P. Bergman) is based on the work by Rybicki & Hummer (1991) and Rybicki & Hummer (1992) and the bench-marking has been reported by Maercker et al. (2008). It takes line overlap explicitly into account and, in the following, we use this program to model the multitude of hfs lines of N₂H⁺. The geometry is spherically symmetric, which needs to be kept in mind when comparing model line profiles with observed ones. Dust continua are included, but switched out in the simulated on-off observations presented below (zero baseline).

The radiative transfer in the hfs transitions with line overlap was computed using the energy levels that we derived from the frequencies presented by Pagani et al. (2009). The collisional de-excitation rate coefficients were adopted from Lique et al. (2015), who have calculated the rates for collisions with H₂ for the temperature range 5 to 70 K and for the lowest hfs levels. However, F. Lique kindly made available to us their entire results for upper energy levels up to J_up = 7 prior to publication. The computation of the inverse rates, that is, the collisional excitation rate coefficients, assumed thermal balance.

The region was devided into three parts: the north, centred on N 2, the middle, centred on SM 1, and the southeast, centred on N 6, where the N-nomenclature refers to the source designations by Di Francesco et al. (2004). However, only pixels with ≥ 3σ integrated intensities in both the N₂H⁺ (3−2) and (6−5) lines were included. These regions, assumed to be spherical, are about 75′′ in size, hence neighbouring regions are overlapping and radii are 0.02 pc (6.75 × 10¹⁶ cm). We kept the turbulent velocity constant throughout the core, that is, υ_turb = 0.25 km s^-1. After a number of initial trials spanning over a much larger parameter space, the ALI code was run on a finer grid around the most probable values, defined by 4.0 ≤ log n (H₂) ≤ 7.0 in steps of 0.2, −11.0 ≤ log X(N₂H⁺) ≤−8.0 by 0.2 and 5 ≤ T_kin ≤ 50 K by increments of 1 K, generating 11776 models. Among these, 2816 solutions were found that fit the observed intensity ratios within the observational errors. Using the individual line profiles, that is, both the shape and the absolute intensity, for both lines, N₂H⁺ (3−2) and (6−5) from the same position (pixel), the best-fit ALI solution was selected for that pixel. As best-fit-estimator, we used sums of the residuals, where the higher frequency spectra were weighted more strongly to compensate for their overall lower intensities. In other words, we minimized , resulting in a map of best-fit model profiles that, together with the observed ones, are displayed in Fig. 4.

In Fig. 5, the spatial distributions of the best-fit model parameters N(H₂), T_kin, n(H₂) and X(N₂H⁺) are shown. Values for the 10′′ pixels span two orders of magnitude in column density, N(H₂). For the star ρ Oph, Bohlin et al. (1978) determined the anomalous relation N(H I + H₂) /E(B−V) = 15.4 × 10²¹ atoms cm^-2 mag^-1. Since R_V = A_V/E(B−V) = 5.5, the gas-dust relation for ρ Oph reads N(H₂) = 1.4 × 10²¹ A_V cm^-2. Using this calibration also for the core, the ALI results over 30′′ scales would imply the surprisingly low visual extinction, A_V, through ρ Oph A of about 100 mag or lower. Average column densities of a few times 10²² cm^-2 are in accord with previous determinations.

Kinetic gas temperatures are about 10−13 K in the central regions, but 17−19 K in the north, and16−19 K in the southeast. Overall, these temperatures are lower than what has been determined before (e.g., Bergman et al. 2011, and references therein). The relative abundance, X(N₂H⁺) = N(N₂H⁺)/N(H₂) appears to peak in the central part, where log X(N₂H⁺) =−9.0 ± 0.2, whereas X(N₂H⁺) is generally lower than 10^-10 in the outer regions. There, derived overall gas volume densities, n(H₂), are around 10⁶ cm^-3, which would be in accord with what might be required by some models of the dust opacity (see below, Sect. 4.4.2). Toward SM 1, and on much smaller scales, Friesen et al. (2014) estimated densities in excess of 10⁹ cm^-3.

As expected, optical depths in the centre of the (3−2) main line are high, explaining the low contrast between main and satellite components (Fig. 5). In contrast, the (6−5) lines are optically thin, which is supported by their observed Gaussian line shapes. The figure also shows the distribution of the radial velocity residuals, Δv = υ₀−υ_LSR, where υ₀ is the observed central velocity and υ_LSR = +3.5 km s^-1. Locally, velocity shifts are lower than 0.2 km s^-1, which is much lower than the assumed turbulent velocity of the gas (0.6 km s^-1). On much larger scales, > 40′, Liseau et al. (1995) derived a NW-SE radial velocity gradient of about 1 km s^-1 pc^-1. For the limited extent of the region studied here, this has only a negligible effect on the line transfer.

4.2. Dust temperature, opacity, and mass

4.2.1. Point-source fitting of the spectral energy distribution (SED)

The fluxes of several point sources in the field were extracted and corrected for the underlying cloud emission. Their SEDs were determined assuming modified blackbody emission, where the adjustable parameter β is the power-law exponent of the frequency dependence of the flux in the Rayleigh-Jeans regime (see Appendix A). The fitting results are presented in Table 3 and shown in Fig. 7. These suggest a temporal evolution of the parameter β, in the sense that the FIR/submm SEDs become increasingly flatter with time. At early times, Class −1 and Class 0, the FIR/submm emission is entirely dominated by the envelope, whereas at the later stages, the envelope emission has essentially vanished and is dominated by the disk or the remnants of it (Class III). A mixture of these two dust-emission components is exhibited by the objects in between (Class I and II). The value of β is presumably dominated by the size distribution of the dust grains (e.g., Miyake & Nakagawa 1993; Kruegel & Siebenmorgen 1994).

To examine the β-dependence on the size a of (spherical) grains, we performed Mie-calculations (Min et al. 2003) using laboratory measurements for a few typical grain materials (Fig. 8). This plot is qualitatively similar to Fig. 4 of Testi et al. (2014) for young disks, but here, β has been calculated for the wavelength region 200 μm to 1.3 mm instead, and it falls into the interval 0.0 ≤ β ≤ 2.5. The power-law exponent p has been varied within the range −5.0 ≤ p ≤ −2.0, and where dn(a) ∝ a^pda. As expected, the β-dependence on chemical composition decreases with increasing p, that is, grain size, and these larger grains are found around the Class II and III objects, viz. LFAM 3 and GSS 32, respectively, with a ≳ 100 μm. In contrast, the very small grains are found predominantly around the probably youngest object, that is, the starless clump SM 1 (Class−1). Around the Class 0 to I sources VLA 1623, GSS 30 IRS1 and LFAM 1, sizes would fall roughly within 0.1 ≤ a ≤ 10μm.

Therefore, the observed time dependence suggests that grain growth mainly occurs in the protoplanetary disks, with only little enrichment of large grains in the parental cloud. However, the precise dating of these events is difficult, but they probably fall within the range of some thousand years (starless cores, Class -1) to some megayears (near the ZAMS, Class III). If our observational results apply to currently adopted pre-main-sequence timescales (e.g., Evans et al. 2009), then β(t) can be fit by a (broken) power law, β ~ t^{− γ}, where γ(<10⁵ yr) ~ 0.04 and γ(≥10⁵ yr) ~ 0.5 to 0.6, respectively. If β scales inversely with grain size/mass, this may be consistent with stochastic processes in a turbulent medium (see also Blum 2004). Apparently, very little happens up to some 10⁵ years. Thereafter, grain growth seems to occur at a faster rate.

Needless to say, it would be very valuable to extend this type of analysis to other star-forming regions, both for better statistics, but above all also for the comparison among different environments. β-evolution around young sources has recently been found elsewhere, but generally only for a single object at a time or, in one case, for two sources (Testi et al. 2014, and references therein). However, the very compactness (0.05 pc ~ 10⁴ AU) of an entire star-forming region in combination with the multitude of objects of different type as discussed here, makes the present study of ρ Oph A unique.

As for SM 1, these results refer to source scales of some 10³ AU. On smaller scales, that is, a few 10² AU, SM 1 appears to have fragmented into smaller pieces (Nakamura et al. 2012). This may be expected on the basis of current theoretical models of star formation (Lomax et al. 2014, and references therein). However, using ALMA at higher resolution, Friesen et al. (2014) found no evidence for substructures down to scales of some 10 AU.

4.2.2. Extended dust emission

The methods for estimating dust emissions from point and extended sources differ only slightly (e.g., Andre et al. 1993), and in Appendix A, we present the analysis of our submm data assuming conventional methods, again for the purpose of comparison with other works. We did not, however, as is often done, assume a single dust temperature nor a single power law exponent β of the frequency dependence of the dust opacity, nor a single normalization value at some fiducial frequency. Instead, we intend to exploit observations at, at least, three frequencies to determine the spatial distribution of the dust temperature, the dust optical depth, and the absolute grain mass opacity. However, from the observed spectral index α = Δlog I_ν/ Δlog ν for the wavebands 350 μm–870 μm and 870 μm–1.3 mm, respectively, it is apparent that the APEX data are internally consistent, but appeared incompatible with the IRAM data (Fig. 3). Whereas α generally is > 2 for the Saboca-Laboca maps, as one would expect for thermal emission in the Rayleigh-Jeans (RJ) regime, the index is < 2 everywhere for the Laboca-IRAM data. We have to conclude that the relative calibration of the LABOCA and SABOCA maps is good, but that the LABOCA and IRAM data suffer from inter-calibration uncertainties, and the latter were therefore not used in our analysis (see Appendix A).

With I_ν obtained from observation, the right-hand side of Eq. (A.2) contains two unknowns. For the interpretation of observations of molecular clouds, it has been common to assume ad hoc values for κ_ν and T_d. In addition, one also assumes that the run of κ with ν is given by a power law, with an exponent β, that is, κ_ν ∝ ν^β. For thermal emission in the RJ-regime, β would be related to the spectral index α, given by β = α−2, and for spherical Mie particles, β takes values between 0 and 2 (e.g., Emerson 1988), but values lower than 0 and greater than 2 are possible when the proper physics of non-perfect crystalline material are included (Meny et al. 2007). The use of the β-relation for mass estimates requires the normalization at some fiducial frequency (cf. Hildebrand 1983). In many cases, the normalization value of the opacity incorporates an assumed g2d of 100 (Beckwith et al. 2000, and references therein).

The normalization is commonly taken from published works on dust opacity computations. However, by using these theoretical results directly, one circumvents the necessity of specifying β-values, which can change from one frequency band to another. Instead, one would wish to determine how these might vary spatially across a cloud, which might tell us something about the cloud dust properties. Potential problems include noise fluctuations and calibration uncertainties, and gradients in density and/or temperature along the line of sight (LOS) that may lead to non-convergence/inconsistent κ-values. Similar would apply to derived β-values, when one assumes that κ(ν) is given by a single power law of index β (e.g., Juvela et al. 2013). The accuracy of our results below is limited to twice the statistical error of the intensity I_ν at each frequency ν. However, the largest uncertainty comes from the absolute value of the opacity.

As an example, we chose the work by Ossenkopf & Henning (1994), which provides dust mass opacities for conditions thought to prevail in dense molecular clouds. Processes of coagulation are considered and the steady state is reached typically after some 10⁵ years, which is most likely younger than the ages of molecular clouds. The initial state consists of bare grains, distributed in size as in the diffuse ISM (viz. “MRN-distribution”, Mathis et al. 1977)⁶. With time, these grains grow in size by acquiring ice mantles of various thickness and opacities, between the wavelengths of 1 μm and 1.3 mm, are provided for a range in environmental density, from n(H) = 0 to 10⁸ cm^-3.

Shown below, in Sect. 4.2.4 (Fig. 16), are the results for fixing the opacity curve, that is, No. 22 in Table 4. The distribution of the dust temperature exhibits an average value of 19 K. Although temperatures are around 20 K in most parts of the map, regions with T_d> 25 K exist along a north-south ridge, in addition to a couple of localized shocked regions, near VLA 1623 and GSS 30 (cf. Fig. 1). The dust optical depth at 350 μm, τ_{350 μm}, is much lower than unity, justifying a posteriori the assumption of low optical depth. This is then of course also simultaneously justified for the LABOCA data (maxτ_{870 μm} = 0.005).

A major drawback of the presented analysis is that it ignores the details of the processes that are heating the observed dust. In other words, the energy balance equation is not solved, which neglects the existence of gradients in density and radiation field strength, and hence temperature, in the source. These shortcomings are considered in the physical models described below, which are two 3D radiative transfer models that follow two different approaches. Regarding spatial coverage, these models are complementary, meaning that model I provides a larger, global view and model II focuses on the central regions of ρ Oph A.

4.2.3. 3D dust modelling: Model I

Model I uses the code RADMC-3D⁷ developed by Dullemond and benchmarked by Pinte et al. (2009). The model consists of a low-density widespread cloud medium, the high-density core ρ Oph A with the dust clump SM 1, the two point sources VLA 1623 and GSS 30, and the two early-type stars S 1 (B 4) and HD 147889 (B 2.5), each with its own PDR bubble. These are the main energy sources in the model, providing the heating of the dust. The stellar SEDs are based on Kurucz ATLAS 9 model atmospheres with the parameters of Table 5.

In detail, the model assumes a high-density core of radius 0.045 pc that is embedded in a low-density cloud. The early-type stars S 1 and HD 147889 generate two PDRs, viz. west and east of the dense core, with respective radii of 0.06 pc and 0.5 pc. The bright B2 star HD 147889 is situated more than half a parsec away to the southwest and behind the cloud (Liseau et al. 1999). Together, these stellar radiation sources are heating the dust of ρ Oph A from the outside. By comparison, internal heating sources are essentially negligible.

The exploited grain opacities are shown in Fig. 9, and in ρ Oph A, individual sources contribute to the total modelled dust mass of 10.4 × 10^-2 M_⊙ as follows: the core with 7, SM 1 with 1.7, VLA 1623 with 0.9, and GSS 30 with 0.8 in units of 10^-2 M_⊙. A match to the PDR around HD 147889 assumes 0.04 M_⊙ of dust (see Fig. 10). Based on C¹⁸O and other molecules (Paper II, in prep.), the density profile is given by ρ ∝ r^-1.5 for ρ Oph A globally and by ρ ∝ r^-2 for SM 1 locally.

In the model, the stellar object S 1 (Chini 1981, B5)⁸ is situated off-centre in a low-density bubble of radius 0.06 pc that is filled with 7 × 10^-4 M_⊙ of 0.1 μm sized carbon particles⁹ (Dorschner et al. 1995). These small particles scatter the short-wavelength radiation extremely well, and this is needed to account for the observed intensity distribution in the near- and mid-IR. However, some finite, non-zero filling of the bubble is also required for numerical reasons. For an assumed g2d = 100, the total mass from model I of the entire ρ Oph A core is 23.7 M_⊙ (Fig. 11).

As already mentioned, we used RADMC-3D to solve for the temperature distribution, which is shown along two different cuts on the right-hand side of Fig. 10. Obviously, the grain temperature also depends on its size, and the distributions for the three groups of grains are shown individually. For μm- and mm-size grains, temperatures of the dust amount to 10−30 K near the dense core SM 1. These model results can be used to generate synthetic images of ρ Oph A. These are shown in Fig. 12, together with the observations at six wavelengths, from the mid-infrared to the submm. Evidently, model I reproduces the observed intensity distributions very well, in particular the large-scale features.

4.2.4. 3D dust modelling: Model II

Thermal equilibrium models The modelling follows the procedure outlined by Juvela et al. (2011), which starts by deducing the 350 μm optical depth distribution for a fixed β = 2 from the observed APEX (350 μm and 870 μm) maps. The optical depth map is used to set the column density distribution of the model. However, the LOS density distribution is not directly constrained by the observations. We assumed a Gaussian distribution with a FWHM of 0.05 pc towards the densest part of the core, increasing up to 0.07 pc towards the low-density regions. With these values the LOS extent of the core is similar to its plane-of-sky dimensions. The resulting 3D density distribution was further modified to include a cavity around the S 1 star. The cavity was modelled as a paraboloid with the star S 1 at the focus. Inside the cavity the density was set to 10^-5 times the highest density of the model core. The models represent a box with a size of 0.225 pc at a distance of 120 pc and use 256³ cells, yielding a resolution of approximately 180 AU, or 15. The three-dimensional density distribution and the plane-of-sky surface density distribution are illustrated in the bottom left and right panels of Fig. 13, respectively.

The model was embedded in an interstellar radiation field with intensity three times that of the solar neighbourhood. However, irradiation of the core was dominated by the individual stars, mainly S 1 and HD 147889. The stellar object S 1 was modelled as a 16 000 K, L = 1100 L_⊙ blackbody, while HD 147889 was given a temperature of 20 000 K and luminosity L = 4 500 L_⊙. Because HD 147889 is more than 0.5 pc away from the dense parts of the core compared to less than 0.1 pc for S 1, the latter is the main heating source in most parts of the model core. In addition, the model included two deeply embedded stellar objects with luminosities of approximately 2 L_⊙, but these do not have any significant effect on the radiation field except in their immediate surroundings. The radiative transfer calculations were made using the CRT program (Lunttila & Juvela 2012). The program calculates the dust temperature taking into account the dust self-absorption and heating, and integrates along the LOS to produce maps of dust FIR emission. For comparison with observations, the synthetic maps were convolved with Gaussian kernels corresponding to the APEX beam sizes.

We have calculated grids of models with varying dust properties and density scaling, but keeping the same radiation sources and the three-dimensional density distribution. In model II.1 we adopted the mass extinction from curve No. 22 in Table 4, taken from Ossenkopf & Henning (1994); however, the opacities are given only for wavelengths λ ≥ 1.0μm and no information is given about scattering properties. Therefore, for the shorter wavelength extinction, the model Milky Way A for R_V = 5.5 of Weingartner & Draine (2001) was used, matched to the former at the wavelength of 1.0 μm. The scattering properties were obtained by assuming the same albedo and scattering asymmetry parameter as in the Weingartner & Draine (2001)R_V = 5.5 model. To study the effect of β, in addition to the unmodified β ≈ 1.9 dust, we used models where the extinction curve at wavelengths λ ≥ 100μm was modified by multiplying the Ossenkopf & Henning (1994) absorption cross-sections by (λ/ 100 μm)^a with a = −0.1, 0.1, 0.2, and 0.3, corresponding approximately to β = 1.8, 2.0, 2.1, and 2.2. The extinction curves are shown in Figs. 13a and b. The peak dust surface density was varied in eight steps between 0.00501 g cm^-2 and 0.0130 g cm^-2.

In model II.1, the best fit was obtained with β = 2.0 and a peak dust surface density Σ_dust = 0.010 g cm^-2. Dust models with lower β predict a too high 870 μm to 350 μm brightness ratio. Assuming a gas-to-dust mass ratio of 100, the mass of the main part of the core (defined as having at least 5% of the peak surface density and corresponding roughly to the extent that was used for the 1D analysis in Sect. 4.2.2) is 7.3 M_⊙. The mass included in the whole simulated region is approximately 20 percent higher. Although the best models fit the observed 350 μm and 870 μm intensity within approximately 20 percent, none of the models is successful in reproducing the 100 μm and 160 μm Herschel-PACS observations. The calculated surface brightness is at least a factor of two higher than observed. In particular, it is difficult to simultaneously produce the relatively high brightness at 350 μm and the low observed intensity in the PACS bands. Figure 14 shows the comparison between the simulated maps from the best-fit model II.1 and the Herschel-PACS 100 μm and 160 μm, and the SABOCA 350 μm, and LABOCA 870 μm maps. To study the effects of numerical resolution on the results, we increased the resolution by a factor of four, that is, to a cell size of 0.0002 pc = 45 AU. The differences between lower and higher resolution were lower than 10% everywhere and lower than 2% where the processes were active. In the following, we therefore ran the models with the lower resolution to keep the computations tractable.

In model II.2 we have chosen BARE-GR-S of Zubko et al. (2004, Table 7) as our basic dust opacities¹⁰. As with model II.1, we ran simulations with modified dust models with different values of β. The modifications of the extinction law were similar to model II.1, except that the changes were applied at wavelengths λ> 20μm. In addition to the standard β ≈ 2.0 model we used β = 1.5, 1.8, 2.1, and 2.2. The density scaling was varied with nine steps corresponding to peak dust surface densities (at 15 scale) between 0.0197 g cm^-2 and 0.0985 g cm^-2. As the Zubko et al. (2004) dust models have much lower sub-mm opacities than the dust used in model II.1, in both cases the models span a similar range of 350 μm optical depths. The best-fit II.2 model has a peak dust surface density Σ_dust = 0.0498 g cm^-2 and β = 2.1. Modified dust with a steeper sub-mm extinction law matches the 870 μm to 350 μm brightness ratio better than the original dust model in this case as well. Using the same A_V-to-gas mass ratio as in the BARE-GR-S model of Zubko et al. (2004) for our modified β = 2.1 dust, the total mass of the best-fit model is 59.1 M_⊙, where a g2d = 161 has been assumed.

In all models, the amount of dust in the cavity is very low, resulting in a visual extinction between the exciting source S 1 and the inner edge of the cavity of A_V = 3 × 10^-3 mag. The extinction toward the observer is less well constrained, but model II.2 gives A_V≲ 15 mag, with similar values for the others. This value compares very well with the A_V = 12.5 mag determined by Ward-Thompson et al. (1989) from the H₂O-ice absorption feature at 3.1 μm, and the range A_V = 11 to 13 mag, for R_V = 5.5, based on the X-ray analysis by Vuong et al. (2003). Non-equilibrium models Models II.1 and II.2 assume that the dust grains are in thermal equilibrium with the local radiation field. However, absorption of a single high-energy photon can heat a very small dust grain to a high temperature. Absorbed energy is therefore reradiated more in the infrared than in the sub-mm (see, e.g., Draine 2003). To examine whether an accurate calculation of the emission from very small grains changes the results, model II.3 uses the dust model Milky Way A for R_V = 5.5 of Weingartner & Draine (2001) but does not assume thermal equilibrium. Because these calculations are much more time-consuming, only three different density scalings were used with this dust model. At wavelengths λ> 200μm the calculated emission from the best-fit model is almost identical to the best-fit model II.2. At 100 μm and 160 μm, model II.3 produces slightly lower fluxes as shown in Fig. 15, but the discrepancy between the model and Herschel-PACS observations remains large. The best-fit model has a peak dust surface density of Σ_dust = 0.0491 g cm^-2 and, using a gas-to-dust mass ratio of 105.1, a mass of 38.0 M_⊙.

The LOS mass-averaged dust temperature, T_dust, and the 350 μm optical depth, τ_{350 μm}, from the best-fit 1D model (see Sect. 4.2.2) and the three 3D models II.1, II.2, and II.3 are shown in Fig. 16. The results from all three 3D models are similar. Their LOS optical depths agree within 20 percent and the calculated dust temperatures only differ by approximately 1 to 2 K. The large differences in calculated masses are caused mostly by the widely varying sub-mm opacities and, to a smaller degree, by the gas-to-dust mass ratios of the dust models. The LOS average temperature through the dense core is between 18 and 20 K, but along the borders of the low-density cavity around the S 1 star the temperature rises to 25–30 K. Results from the 1D model show much larger variations in temperature, although the average value is similar. The 350 μm optical depths from the 1D analysis are approximately 30% lower.

4.3. The g2d in ρOphA

Table 6 compares the results of different approaches to the dust modelling. Derived absolute masses can vary quite substantially, from 5 to 40 M_⊙. However, opacity-weighted, normalized masses agree to within 10% of each other, with an average value of 1.1 over an interval of 0.4 to 1.4.

To determine the absolute gas-to-dust mass ratio, g2d, is therefore not a trivial task. For example, if we choose the result of the non-equilibrium model 2.III, the map average is only slightly sub-canonical, g2d = 88 (Fig. 17). This value would directly agree with that derived from X-ray work by Vuong et al. (2003), who found g2d = 95−80 within the range R_V = 4−6¹¹.

However, on scales of the individual sub-regions, this number would not be truly representative: over tens of seconds of arc, g2d takes median values of around N 2 in the north, of in the central regions of SM 1, and of in the southeast, near N 6. Taken literally, local g2d values in ρ Oph A imply fractions of the dust mass from lower than 0.5% to more than 10%. In particular, the deviating, low values around the core SM 1 were already remarked upon by Bergman et al. (2011).

To bring the central g2d into line with the widely accepted value of 100 would require the downward adjustment of the derived gas mass by a factor of ten. This is not a realistic option on the basis of the ALI modelling, however, and must be discarded. Similarly, reducing the dust mass by an order of magnitude would result in g2d of several thousand in the north and southeast, an option that is not convincing either. On the other hand, an uneven projected distribution of gas and dust may be a solution. In Fig. 18, maps of gas and dust at 7′′ resolution are displayed, from which it becomes apparent that the continuum radiation peaks at positions offset from those of the gas. A low g2d would naturally result.

4.4. Distribution of gas and dust in ρ Oph A

The C¹⁸O (3−2) line and the LABOCA continuum band are close in frequency, at about 850 μm wavelength, and the maps of both have been observed with the same telescope (APEX), meaning that these data share the same angular resolution of ~20′′. Our SABOCA continuum map at 350 μm has also been obtained with APEX, albeit at the higher resolution of 7′′. Both APEX continuum maps share the same offset centre and internally agree excellently well (Sect. 4.2). The nominal offset of the C¹⁸O (3−2) map (Liseau et al. 2010) is at (−451, + 30), which was applied to the line data to align with the continuum maps. As already mentioned, the pointing accuracy of the APEX telescope is 2′′ rms.

We applied the program IMAGINE, which is based on an in-house developed algorithm (Rydbeck 2008), to the spatially oversampled C¹⁸O (3−2) data. The number of noise channels in the line spectra, in each positional pixel, is large, which allows an accurate description of the statistics of the noise distribution. In addition, the spectra have a high signal-to-noise ratio. These properties of the data are used by IMAGINE to deconvolve high-quality data, increasing the resolution by up to a factor of three. The deconvolved C¹⁸O (3−2) map of the integrated line intensity for the gas has then effectively the same resolution as the 350 μm continuum data for the dust (Fig. 18). It is gratifying to notice that the point source VLA 1623, which is detected in both data sets, is at the same position in both maps. Furthermore, as a final sanity check, we degraded the high-resolution result with the APEX beam at 850 μm, and this procedure did indeed recover the observed (input) map. We are therefore confident that any differences in the intensity distributions and on spatial scales exceeding 3σ ( = 6′′) of the pointing accuracy are real.

This gas-dust dichotomy is not immediately apparent in the presented maps involving N₂H⁺, since these have a much lower resolution, governed by the (J = 6−5) observations with the roughly 40′′-beam of Herschel. However, the (3−2) data have been obtained with APEX and for the same coordinates as the LA-/SABOCA maps. Using IMAGINE for the N₂H⁺ (3−2) map then yields a similar resolution as the SABOCA image and the sharpened C¹⁸O (3−2) map. The result is shown in the right panel of Fig. 18. The clear offset, by 12′′, between the continuum source VLA 1623 and an intensity peak in N₂H⁺is at first astonishing. Similar mismatches are also found elsewhere in these images, and they have also been noticed by others (Di Francesco et al. 2004; Friesen et al. 2014). We therefore conclude that the spatial non-coincidence between moleculer gas and dust is real. However, since the N₂H⁺ emission ridge is much narrower than that in C¹⁸O, a depression in the distribution is not evidenced by the present image. This may require data at higher resolution.

From Fig. 18, it is immediately evident that the dust continuum and the molecular line emission have different spatial distributions in the sky. Maximum line emission in C¹⁸O (3−2) is observed in a ridge closest to the stellar heating source S 1, whereas the continuum is mostly pronounced offset to the west, where all dust models, both 1D and 3D, indicate dust temperatures T_dust ≥ 20 K. In contrast, the spectral line analysis points toward gas temperatures that are clearly much lower than twenty Kelvin (see also, e.g., Stamatellos et al. 2007). This might effectively have created the “hole” of the C¹⁸O emission, which might be due to depletion of the CO gas onto the dust. At T_kin< 25 K, carbon monoxide is mostly frozen out (Sandford & Allamandola 1993; Öberg et al. 2005), whereas diazenylium is possibly much less affected because it mainly remains in the gas phase. It is not clear, however, how N₂H⁺ would react at temperatures even below 15 K (Öberg et al. 2005). A selective freeze-out scenario would not explain the low g2d that is based on N₂H⁺ observations. On the other hand, if all molecules are frozen out, there would be no contradiction¹².

5. Conclusions

In the FIR/submm, observed point sources exhibit a dependence of the dust parameter β on evolutionary status (SED class), β = β(t). Values decrease monotonically with age. For the pre-stellar core SM 1, the continuum emission is entirely dominated by the thick envelope and β = 2, whereas in pre-main-sequence stellar objects of Class III, the envelope is lost and the emission originates only from the circumstellar disk, β ~ 0. Objects of intermediate type display intermediate β-values. If not merely a matter of contrast, this may imply that grain growth primarily occurs in the protoplanetary disks, not in the natal cloud cores, since low values of β are ascribed to large grains.

ρ Oph A is exposed to uneven irradiation by stellar sources, heating the core primarily from the outside. We examined different dust emission models, including both 1D and 3D geometries, without and with thermal balance and with non-equilibrium conditions. Mass estimates span a factor of seven, but do agree to within a factor of two when proper weighting by the absolute mass opacity κ_ν is applied.

To model the observed gas emission in N₂H⁺ lines, best-fit solutions were obtained by combing a multi-dimensional parameter space. These provided excellent fits to the spatial distribution of the line intensity and the line shape. Combining the results for the gas with those for the dust yielded a global gas-to-dust mass ratio that is only slightly sub-canonical, g2d = 88. However, on local scales, this hardly applies anywhere, with ratios much higher than 10² in the north and south and significantly lower than that in the central regions. There, gas and dust appear to be spatially segregated. At high angular resolution, a hole of C¹⁸O (3-2) emission is accompanied by maximum emission in the dust continuum. At lower resolution, this may result in a misleading g2d, and perhaps indicate that N₂H⁺ is frozen out as well.

---

Acknowledgments

We thank the anonymous referee for the very thoughtful report that led to a number of clarifications in the manuscript. François Lique kindly made available to us the collision rate data for N₂H⁺ prior to publication, which we hereby gratefully acknowledge. R.L. enjoyed interesting discussions with John H. Black. The help by the APEX staff with our ground based observations is very much appreciated. We thank the Swedish National Space Board (SNSB) for the continued support of our Herschel projects. The computations were performed on resources at Chalmers Centre for Computational Science and Engineering (C3SE) provided by the Swedish National Infrastructure for Computing (SNIC).

References

Appendix A: Estimation of dust masses in homogenous and isothermal media

The modification we made to the commonly adopted procedure is based on the fact that we do not a priori assume only a single dust temperature and grain opacity. Rather, these quantities are parameters to be fit iteratively by comparison with submillimetre observations. The value of the optical depth at any frequency is obtained and tested a posteriori. The specific intensity I_ν, in erg cm^-2 s^-1 Hz^-1 sr^-1, along the LOS of optically thin thermal emission from dust particles can be written as (A.1)where κ_ν is the dust mass opacity, in cm² g^-1, at the frequency ν, in Hz, n_d is the volume density of dust particles, in cm^-3, m_d is the mass of a dust grain, in g, and B_ν(T_d) is the Planck function, in erg cm^-2 s^-1 Hz^-1 sr^-1, at ν and for the dust temperature T_d, in Kelvin.

For constant κ_ν, m_d and T_d along a ray, the integral is only over the volume density n_d, meaning that it is the column density of dust grains, . The surface density of dust particles is Σ_d = m_dN_d. Identification in Eq. (A.1) yields (A.2)From Eq. (A.2) it is apparent that the optical depth through the dust is given by (A.3)The assumption that τ_ν ≪ 1 will have to be justified a posteriori. Applying Eq. (A.2) to observations I_ν1 at a given frequency ν₁, together with κ_ν1 from a particular opacity curve (j = 1, 2, 3,....), yields a multitude (i = 1, 2, 3,...) of solutions for both Σ_d and T_d, that is, (A.4)This multitude is limited by the range of examined temperatures, viz. T_d = 4 to 2000 K, where the first value refers to the average equilibrium temperature of grains irradiated by the interstellar radiation field in the solar neighbourhood¹³ and the second to approximate dust evaporation temperatures. Similarly, the number of discretized temperature values, T_d(i,j), that need to be computed is ultimately limited by the measuring accuracy of I_ν.

This multitude can be further reduced, in fact to unity, since, for any given opacity curve, only one value of T_d (and Σ_d) will simultaneously satisfy Eq. (A.2) for data at another frequency, for example, I_ν2, viz. (A.5)where we let ΔI_ν2(j) → 0. Equation (A.5) yields for each opacity curve, designated by j, a unique solution for T_d(j) and Σ_d(j). It only remains to select the correct κ-curve. Minimization of the data for a third frequency, I_ν3, could in principle be used to achieve this. Hence (A.6)Minimizing Δκ_ν then provides the unique solution for the temperature T_d, the surface density Σ_d and the opacity law of the dust, κ = κ(ν). Finally, the assumption of low optical depth at the different frequencies of this solution is then verified a posteriori (Eq. (A.3)).

The column density of the gas, N_g = N(H), needs to be determined from other observations. This, finally, leads to the gas-to-dust mass ratio from

(A.7)where μ = 2.4 is the mean molecular weight. The derived dust opacity curve has been computed for a particular volume density of the gas, n^d(H) and this ought to be consistent with that determined from the spectral line fitting, n^g(H). The spatial averaging is done over each pixel, so that the total dust mass, M_d, is then obtained from the summation of the surface density per pixel over all pixels.

All Tables

All Figures