© ESO, 2016

1. Introduction

Apart from in situ measurements in the solar system, all information gained in astrophysics is deduced from the analysis and interpretation of radiation received with ground or space-based telescopes. Gas and dust, in between the source and the telescope, influences the radiation. This means that analyzing the radiation received from an astrophysical object not only provides information about the source, but also about the medium in between the object and the observer. Radiative transfer is thus one of the most fundamental phenomena in astrophysics.

As summarized by Steinacker et al. (2013), three-dimensional dust radiative transfer calculations are essential for making progress in many fields of astronomy. Dust grains modify the radiation field in many objects, such as protoplanetary disks, evolved stars, reflection nebulae, supernova remnants, molecular clouds, the interstellar medium, galaxies, galactic nuclei, and the high-redshift Universe.

In this paper, we focus on the modeling of the dust and free-free continuum emission of the high-mass star forming molecular cloud Sagittarius B2 (Sgr B2) by applying detailed 3D radiative transfer modeling.

Sgr B2 is one of the most massive molecular clouds in the Galaxy. It is located at a distance of 8.34 ± 0.16 kpc (Reid et al. 2014)¹ and has a projected distance of 107 pc (or 43.4′ ) from Sgr A^∗, the compact radio source that is associated with the supermassive black hole located at the Galactic center. Hüttemeister et al. (1993) distinguish three different parts in Sgr B2: (i) a low density envelope; (ii) a moderate density region extended around (iii) local hotspots, which are the most compact, densest molecular regions (see Fig. 1). The envelope measures 38 pc (or 15.4′ ) in diameter (Scoville et al. 1975, corrected for distance) and has a gas mass of 7 × 10⁶ M_⊙ (Goldsmith et al. 1990). The average H₂ density n(H₂) ~ 10³ cm^-3, and H₂ column density N(H₂) ~ 10²³ cm^-2, are relatively low compared to the central part of the region. The moderate density region extends over 2.5 pc × 5.0 pc around the local hotspots. Its density and H₂ column density are higher, n(H₂) ~ 10⁵ cm^-3 and N(H₂) ~ 10²⁴ cm^-2. Embedded in this intermediate region are the local hotspots. At least three of them are sites of active star formation (Gordon et al. 1993). These three sources are historically named according to their relative location in an equatorial coordinate system: Sgr B2(N)(orth), Sgr B2(M)(ain), and Sgr B2(S)(outh). They are positioned along a north-south line. In projection, Sgr B2(M) is located 1.9 pc (or 45′′ ) south of Sgr B2(N) and Sgr B2(S) is located 1.7 pc (or 42.3′′ ) south of Sgr B2(M), see Fig. 1. These cores have sizes of ~0.5 pc, H₂ densities of ~10⁷ cm^-3, and column densities of 10²⁵ cm^-2 (Hüttemeister et al. 1995). Different works, e.g., Goldsmith et al. (1990), Etxaluze et al. (2013), have derived the gas mass of the two cores Sgr B2(N) and SgrB2 (M) in the ranges of 6−25 × 10⁴ M_⊙ and 3−23 × 10⁴ M_⊙, respectively. However, they used different radii (0.4−1.6 pc), which makes a comparison of the results impossible.

Sgr B2(N) and Sgr B2(M) are sites of active massive star formation and comprise a plethora of (ultra-compact) Hii regions, X-ray sources associated with Hii regions, X-ray sources with no radio or IR counterparts (Takagi et al. 2002), dense cores, embedded protostars, and molecular masers (Goicoechea et al. 2004). More than 70 Hii regions have been detected in the whole Sgr B2 cloud complex (see Fig. A.1; Mehringer et al. 1993; Gaume et al. 1995; De Pree et al. 1998).

In the observed dust continuum maps, we see an extension of the cloud complex to the north-east. Following the historical naming scheme, we name this region Sgr B2(NE) throughout this paper. This component is located in projection 2.3 pc (or 56.8′′ ) north-east of Sgr B2(N). Another extension is visible towards the south of the cloud. We will name this extension Sgr B2(deepS)outh throughout this paper. This component is located in projection 3.9 pc (or 95.6′′ ) south of Sgr B2(S), see Fig. 1.

The modeling of the continuum emission of the Sgr B2 complex presented here provides us with the three-dimensional model of the structure (density distribution) of this region. For this multiwavelength, multiscale data is crucial to properly constrain the structure. In the next step this will enable us to model the line shapes (kinematics) and thus constrain molecular gas properties such as the gas velocity field, molecular abundances, etc.

The paper is organized as follows: in Sect. 2 we introduce the observational dataset used throughout the paper. The modeling approach is presented in detail in Sect. 3. This is followed by the application of the modeling approach to Sgr B2 and the discussion of the results in Sect. 4. Finally the paper is concluded in Sect. 5. In Sect. A, we present additional figures and in Sect. B we tabulate the setup parameters of all models. In Sect. C we derive physical properties of Hii regions.

2. Observations and data reduction

Multiwavelength, multiscale data is crucial to properly constrain the structure of Sgr B2. Towards the hot cores Sgr B2(N) and Sgr B2(M), the Herschel/HIFI spectral surveys provide the continuum information from the sub-mm up to the far-infrared regime. High-resolution interferometric maps towards both hot cores obtained with the Submillimeter Array (SMA) and the Very Large Array (VLA) provide the necessary spatial resolution on small scales. To cover the large-scale structure, we use dust continuum maps obtained within the surveys ATLASGAL and HiGAL, described in detail below. A summary of the data is presented in Table 1. Figure A.6 provides an overview of the employed datasets and the spatial scales they cover.

2.1. Herschel/HIFI

The Herschel/HIFI guaranteed time key project HEXOS (Herschel/HIFI observations of EXtraOrdinary Sources; Bergin et al. 2010) includes full line surveys of Sgr B2(N) towards α_J2000 = 17^h47^m19.88^s, δ_J2000 = −28°22′18.4″ and Sgr B2(M) towards α_J2000 = 17^h47^m20.35^s, δ_J2000 = −28°23′03.0″, covering the frequency ranges of 480−1250 GHz and 1410−1910 GHz. The corresponding half-power beam widths are 44.9−17.2′′ and 15.3−11.3′′ , respectively.

The spectral scans have been calibrated with HIPE version 10.0 (Roelfsema et al. 2012). The resulting double-sideband (DSB) spectra were reduced with the GILDAS CLASS² package. Basic data reduction steps included removal of spurious features or otherwise unusable parts of the spectra. The continuum emission was subtracted from the DSB scans by mostly zero-th, first- or more rarely second-degree polynomial fitting. The continuum-subtracted DSB data were deconvolved (sideband separation through pure χ² minimization; Comito & Schilke 2002) to provide an equivalent single-sideband (SSB) spectrum for each HIFI band.

Since a single full HIFI line survey is actually made up of 14 independent line surveys with seven different local oscillators (HIFI LO bands 1a through 7b), inconsistencies in the continuum level between HIFI bands were expected and indeed observed. A linear least squares fit of the subtracted continuum values as a function of local oscillator (LO) frequency provided a reliable – because unaffected by spectral features – parametrization of the continuum variation across each HIFI band, which was then folded back into each continuum-subtracted SSB spectrum.

Finally, the overall continuum was rendered self-consistent in two steps: the first adjustment consisted of an additive factor for each band, to achieve a zero-continuum level for the observed saturated absorption features. This is based on the absorption of molecules with a high dipole moment. These molecules are located in foreground, i.e. low density, clouds along the line-of-sight towards Sgr B2. Thus they will have a very low excitation temperature. The second adjustment required a multiplicative factor, in order for the continuum values in overlap regions between bands to be consistent with each other.

For Sgr B2(N), the additive factors range between −1.08 and 0.48 K, with a median of 0.05 K; the multiplicative factors range between 0.82 and 1.42, with a median of 1.00. For Sgr B2(M), the additive factors range between −0.87 and +0.34 K, with a median of 0.28 K; the multiplicative factors range between 0.86 and 1.42, with a median of 1.03.

2.2. Submillimeter Array (SMA)

Sgr B2 has been observed with the SMA in the compact and very extended configurations. The observations were carried out on June 11, 2010 using seven antennas and on July 11, 2010 using eight antennas. The phase tracking centers were α_J2000 = 17^h47^m19.883^s, δ_J2000 = −28°22′18.4″ for Sgr B2(N) and α_J2000 = 17^h47^m20.158^s, δ_J2000 = −28°23′05.0″ for Sgr B2(M). The data reduction and results are described in detail in Qin et al. (2011). The absolute flux scale is estimated to be accurate to within 20%. Both sources were observed in double-sideband mode, and covered rest frequencies from 342.2 to 346.2 GHz and from 354.2 to 358.2 GHz. The line-free channels have been used to reconstruct the continuum image. We have re-imaged the continuum maps with almost uniform weighting, resulting in a beam of 0.37′′ × 0.22′′, position angle (PA = 17.8° with a rms of ~26 mJy/beam.

2.3. Very Large Array (VLA)

To constrain the location and physical parameters of the Hii regions in Sgr B2, we used several different data sets obtained with the VLA.

2.3.1. Gaume et al. (1995) map at 1.3 cm covering Sgr B2

Gaume et al. (1995) observed Sgr B2 with the VLA in three configurations, BnA, CnB, and DnC, between February 1989 and October 1989. The central pointing position is α_J2000 = 17^h47^m20.166^s, δ_J2000 = −28°23′04.76″. The data set and its calibration is described in detail in Gaume et al. (1995). We used an image of the combined DnCnBnA data, resulting in a beam of 0.27 ′′× 0.23 ′′ (HPBW), PA = 70° , with an rms of ~0.38 mJy/beam.

2.3.2. Rolffs et al. (2011a) maps at 7 mm covering Sgr B2(N) and Sgr B2(M)

Sgr B2 has been observed with the VLA in the BnA hybrid configuration at 40.7669 GHz, corresponding to 7.5 mm (project AR687). The phase tracking centers were α_J2000 = 17^h47^m19.902^s, δ_J2000 = −28°22′17.8″ for Sgr B2(N) and α_J2000 = 17^h47^m20.202^s, δ_J2000= −28°23′05.3″ for Sgr B2(M). Sgr B2(M) was observed on 2009 January 31, and Sgr B2(N) on February 1. The data set and its calibration is described in detail in Rolffs et al. (2011a). The continuum was fitted using line-free channels. We have re-imaged the continuum map with almost uniform weighting, resulting in a beam of 0.15 ′′× 0.10 ′′ (HPBW), PA = 52.8° , with a rms of ~0.9 mJy/beam.

2.3.3. De Pree et al. (1998) map at 7 mm covering Sgr B2(M)

De Pree et al. (1998) observed Sgr B2(M) with the VLA using the A configuration on 1996 December 12, and 1997 January 17. The central pointing position is α_J2000 = 17^h47^m20.115^s, δ_J2000 = −28°23′04.02″. The data set and its calibration is described in detail in De Pree et al. (1998). The image we used has a spatial resolution of 0.049 ′′× 0.079 ′′ (HPBW), PA = 11° , with an rms of ~0.5 mJy/beam.

2.4. APEX/ATLASGAL project

Carried out using the Large APEX BOlometer CAmera (LABOCA; Siringo et al. 2009), the APEX Telescope Large Area Survey of the Galaxy (ATLASGAL; Schuller et al. 2009) covers the full inner Galactic Plane at 870 μm with a resolution of 19.2′′ and an rms below 70 mJy/beam. Each position of the inner Galatic plane has been mapped twice with different scanning directions to avoid striping, using the on-the-fly mapping technique. The pointing accuracy is of the order of 4′′ and the flux calibration uncertainty is lower than 15% . A detailed description of the data reduction is given in Csengeri et al. (2014). This map was recently cross-calibrated using the data of the Planck mission and the large-scale structure that is filtered out during the processing of the LABOCA data has been added back into the map (Csengeri et al. 2016). We use this improved map in our study.

2.5. Herschel/HiGAL project

The Herschel Hi-GAL survey (Molinari et al. 2010) provides photometric mid-IR observations at 70 μm and 160 μm using PACS (Poglitsch et al. 2010), and sub-milimeter observations at 250 μm, 350 μm, and 500 μm using SPIRE (Griffin et al. 2010). Sgr B2 was observed as part of the Field 0 observations covering the Galactic center (OD 481; obsids 1342204102, 1342204103). The observations were carried out in PACS/SPIRE parallel mode with a fast scanning speed of 60 ′′ s^-1. The data reduction is described in detail in Traficante et al. (2011). This dataset was cross-calibrated using the data of the Planck mission on the long-wavelength side and the data from the IRAS mission on the short wavelength side. The angular resolutions at the five wavelengths, listed by increasing wavelength, are 10.4′′ , 13.6′′ , 23.4′′ , 30.3′′ , and 42.5′′ (Traficante et al. 2011). We note that the 160 μm PACS map is saturated towards Sgr B2(M) and N and is thus not used in our study.

3. Modeling procedure

We have implemented a framework calledpandora that follows the flowchart shown in Fig. 2 for the continuum modeling. As part of this framework, we employ the three-dimensional, publicly available radiative transfer program RADMC-3D (Dullemond 2012, version 0.39), Miriad for the post-processing, and MAGIX (Möller et al. 2013) for the optimization of the model input parameters. We describe the model setup and physical parameters in the following subsections. RADMC-3D is written in Fortran90 and allows the user to setup their models in a non-invasive way, i.e. via a separate subroutine. In the following subsections we describe different parts of the setup of our model.

3.1. Coordinate system

We use a positively right-handed cartesian coordinate system. This means that the x-axis points west on the sky, the y-axis points to the north and the z-axis points towards the observer. The origin of this system is located at the center of the model, which is chosen to be Sgr B2(M), i.e. α_J2000 = 17^h47^m20.172^s, δ_J2000 = –28°23′04.58″.

3.2. Grid refinement

We want to recover the large dynamic range in spatial resolution from 38 pc, i.e. the diameter of the envelope (Scoville et al. 1975, corrected for distance) down to 100 au (small scale structure around the hot cores and Hii regions). If we would attempt this with a cartesian grid, i.e. a grid where the elements are unit cubes, this would require 800 billion cells. This is not only computationally unfeasible but also unneccessary when e.g. the density distribution is flat.

We use the adaptive mesh refinement technique (originally presented by Berger & Oliger 1984; Berger & Colella 1989) to locally increase the spatial resolution of our numerical radiative transfer simulation. The initial grid consists of 11³ cells. We employ the RADMC option to include a tree-based AMR method (Khokhlov 1998, and references therein). This means that, on a cell-by-cell basis, parent cells are refined into children cells resulting in a recursive tree structure. The resulting grid provides high resolution where needed based on the refinement criteria.

In our approach, two criteria are used to check if a cell needs to be refined: (i) the presence of a dust core center within the cell and (ii) the dust density distribution within the cell. Since we cover a large dynamical range, it is crucial to make sure that small clumps are recognized by the refinement routine. Thus as long as the cell size is larger than twice the radius of a dust core, this cell is refined if the core center is located within the cell. A cell is also refined into eight sub-cells when (1)

where n_center is the density at the cell center and n_i is the density at the center of the faces or the corners of the cell and ϵ is the allowed maximum density difference, which is set to 10% .

The grid is refined to level 13, i.e. at least one initial cell has been refined 13 times. This results in a minimum cell size of 100 au and in total ~68 million final cells.

3.3. Heating sources

RADMC-3D provides several methods to include luminous sources. We choose the option to manually specify the individual luminous sources and assume that all stars are point sources, i.e. their radius is not taken into account. We use a two-step approach. In the first step we account for observed early-type high-mass stars by including the known Hii regions (Mehringer et al. 1993; Gaume et al. 1995; De Pree et al. 1998, see Table B.1 for a complete list of all parameters). Due to a revised distance of the Galactic center we estimated the corresponding parameters, see Sect. 3.5. We specify their position and calculate the luminosity and temperature from the zero-age-main-sequence (ZAMS) type using Table 5 from Vacca et al. (1996). Here we assume that each Hii region is ionized by a single star. We use the RADMC-3D option to assume simple blackbody spectra for each star, by specifying their blackbody temperature. This accounts for stars down to the B0 spectral type.

In the second step we take later spectral types, i.e. stars which cannot produce Hii regions detectable with current observations, with stellar masses between 0.01 M_⊙ and ~19.0 M_⊙, into account. We do this as follows: based on the gravitational potential, our algorithm randomly determines positions for the new stars. For this, we follow the procedure explained in the appendix of Aarseth et al. (1974). A luminosity, which is randomly drawn from the initial mass function (IMF) of Kroupa (2001), is assigned to each star.

We do this as follows: assuming a spherical symmetric star cluster, we specify the radius of this cluster. This radius is determined using a radial histogram of the distributed Hii regions (see Fig. A.3). From Kroupa’s IMF we then calculate the fraction f_HM of the stars in the high mass regime (i.e. between M₁ and M₂; M₁<M₂) to the mass of all stars (i.e. between M₀ and M₂; M₀<M₁<M₂). The lowest mass boundary M₀ is provided, the mass boundaries M₁ and M₂ are set by the high-mass stars that were distributed within the cluster in the first step described above. The total stellar mass of the cluster M_cluster is then calculated from the ratio of the sum of the stellar mass of all high-mass stars distributed between M₁ and M₂ to the ratio f_HM. We then randomly draw stars between M₀ and M₁ from the IMF until the total stellar mass of the cluster, M_cluster, is reached.

The model can contain several of these star clusters. After the stars are sprinkled separately within each of these clusters, the entire envelope is treated as a cluster itself and the procedure is repeated once more, i.e. necessary stars are sprinkled everywhere except in the previously mentioned star clusters until the calculated total stellar mass of the computational domain is reached.

We then convert the stellar mass M to luminosity L using the parameterized mass-luminosity relation from Griffiths et al. (1988): (2)Griffiths et al. (1988) provide two sets of boundaries. One is non-continuous but has a better χ² value whereas the other set, which we choose to use, provides a continuous behaviour of the mass-to-luminosity relation while having a slightly worse χ² value.

We only include luminosity in terms of ZAMS luminosity. The total luminosity determined in the literature, however, includes contributions from accretion luminosity, as well as from high-mass stars without an Hii region (Hosokawa & Omukai 2009). We thus note that the luminosity calculated with our approach is a lower limit. We determine the corresponding effective temperatures of the stars using the fundamental stellar parameters (spectral type, luminosity, effective temperature) compiled by Straizys & Kuriliene (1981). Since we need to provide an effective temperature for each star, the limitation of the tabulated data sets the lower mass limits of the star clusters to 0.2 M_⊙.

3.4. Dust density distribution

The overall density structure is obtained by the superposition of the density profiles of all dust cores, i.e. in overlap regions, the density simply adds up. In each cell j, the density is determined as (3)where i is the index of the dust cores and N is the number of cores. Following Qin et al. (2011), we use a modified Plummer-like profile to model the SMA dust density cores as well as the large-scale envelopes. For a subset of our models, we find it useful to introduce elongated density structures (see Sect. 4). We define them as follows: (4)where n_c is the central density given in H₂ cm^-3 and is given by the Euclidean norm, including scaling factors (5)where r_x,y,z are the components of r and r_0,x, r_0,y, and r_0,z set the size of the cluster core in each of the three principal axes. A spherical symmetric distribution is obtained by setting r₀ = r_0,x = r_0,y = r_0,z. The density distribution is flat inside the radius r₀ and approaches a power-law with an exponent η at r ≫ r₀.

3.5. H II regions

In our model, we consider Hii regions as Strömgren spheres, i.e. as a fully ionized, spherical regions of uniform electron density with no dust. There are more than 70 Hii regions known in Sgr B2 (Mehringer et al. 1993; Gaume et al. 1995; De Pree et al. 1998).

We obtained the high-resolution interferometric maps presented in Gaume et al. (1995) and De Pree et al. (1998) and performed the following analysis steps. Hii regions come in many different shapes, which makes them hard to fit with Gaussians. So in order to match them as best as possible with spheres, we masked everything below three times the rms in the continuum maps and by eye inspection adjusted circles to enclose the observed Hii regions. We then integrated the enclosed flux and assuming optical thin emission and a homogeneous, non-expanding Hii region, we calculated the number electron density n_e, the emission measure EM, and the number of ionizing photons using the following formulas. Here F_ν is the flux density of the Hii region, T_e is the electron temperature, ν is the frequency, D is the distance to the source, and θ_source is the angular diameter of the Hii region. A derivation of these formulas can be found in the Appendix C.

Deriving the number electron density at a certain frequency from the observed flux density assuming optical thin emission will underestimate the number electron density if the Hii region is actually optically thick. This underestimate will then lead to the intensity of the (optically thick) Hii region being underpredicted by the simulation compared to the observed data at this frequency. To account for this discrepancy, we compare the synthetic intensities obtained from our model setup (see Sect. 3.7) with the observed intensities. From the deviation of the intensity levels, we identify the Hii regions which must be optically thick. We then iteratively adjust the flux density values and recalculate the number electron density n_e, the emission measure EM, and the number of ionizing photons until we obtain a resonably good match between the observed and synthetic intensities. As our main goal is to model the dust emission here we need to derive the free-free contribution from the Hii regions. Hence we did not aim at a very sophisticated model of the Hii regions, which would have to include a different geometry and a density structure.

The Hii regions detected by Mehringer et al. (1993) were included according to the parameters stated therein. However, Mehringer et al. (1993) assumed a source distance of 7.5 kpc. We have thus re-calculated the radii and the parameters stated above for a distance of 8.5 kpc. But we have not corrected the values to account for the optical depth effects. Thus we note that these values are a lower limit. All parameters, including the radii and the precessed coordinates, are listed in Table B.1.

3.6. Dust temperature calculation

The dust temperature is calculated self-consistently by RADMC-3D using the Monte Carlo method of Bjorkman & Wood (2001), with various improvements, such as the continuous absorption method of Lucy (1999). Initially, all cells have a dust temperature equal to zero. To determine the dust temperature, RADMC-3D first identifies all sources of luminosity, i.e. the total amount of energy available. This amount is divided into photon packages, which are separately emitted by the individual stars one after another. As the photon packages move through the grid, they interact with the dust present (scattering, absorption) causing them to change their direction or wavelength. In the case of absorption, the photon package is immediately re-emitted in a different direction with another wavelength according to Bjorkman & Wood (2001). However, the luminosity fraction each photon package represents remains the same. Whenever a photon package enters a cell, it increases the energy of that cell and thus also increases the temperature of the dust. Photon packages never get lost, they can only escape the model through the outer edge of the grid. After the last photon package has escaped the grid, the dust temperature is obtained. In summary the dust temperature of each cell is basically the sum of the energy fractions of each photon package passing through the cell. In total we use 10 million photon packages.

The calculated temperature is an equilibrium dust temperature, since it is assumed that the amount of energy each dust grain acquires and re-radiates stays the same. For most cases, this is presumably a very good approximation, because the heating and cooling timescales for dust grains are typically very short, compared to any time-dependent dynamics of the system (Choudhury et al. 2015). Thus transiently heated small grains are not important within the wavelength range considered here.

Note, we only consider stars as heating sources (see Sect. 3.3). The calculated dust temperature in the outer parts of the envelope is on average 12−15 K. Modified blackbody fitting of SCUBA data (Pierce-Price et al. 2000) and Herschel data Etxaluze et al. (2013) yields average dust temperature values for Sgr B2 (along the line-of-sight) of ~20 K. Apart from line-of-sight effects in the modified blackbody fitting, which do not influence the modeling, part of this discrepancy is most likely due to heating by sources not considered here, e.g. cosmic rays or coupling with the turbulently heated gas.

We furthermore make use of the Modified Random Walk (MRW) method that is implemented in RADMC-3D in the simplified form described in Robitaille (2010). This method prevents a photon package from getting trapped in high-density regions by predicting where the photon will go next. This allows RADMC-3D to make one single large step of the photon package, saving the computation time for the otherwise necessary hundreds or thousands of absorption or scattering events.

3.7. Images and post-processing

We use one dust species throughout the model and do not include scattering events. We use the tabulated dust opacity from Ossenkopf & Henning (1994) for dust without grain mantles and no coagulation, as found to best fit the Sagittarius B2 region by Rolffs et al. (2011b). Including the free-free emission we cover the frequency range from 40 GHz up to 4 THz. Continuum maps at various wavelengths of interest are produced with RADMC-3D.

In order to compare the synthetic maps with real obervations, telescope-dependent post-processing is necessary. For this purpose we employ the data reduction package Miriad (Sault et al. 1995). In case of interferometric observations, we fold the simulated maps with the uv-coverage. The imaging is then performed for both maps (synthetic and observed) with the same imaging parameters. In case of single dish observations, the synthetic maps are convolved with the beam of the telescope at that frequency, which we assumed to be Gaussian. Final maps are then produced for the respective telescopes at various wavelengths in intensity units (Jy/beam).

3.8. Fitting procedure

Three-dimensional modeling intrinsically has many free parameters. For this initial study, we have focused on deriving the density field for a fixed dust setup and, apart from by-hand adjustments to account for optical depths effects, fixed Hii regions. Still, the amount of free parameters is large. We have 20 dust core components and 9 components for the large-scale envelope. For each component we modified, in order of increasing importance, the density exponent, a stellar heating source and accordingly a stellar cluster, the position along the line of sight, the displacement in right ascension and declination, the radius and the central density. In total this amounts to approximately 140 free parameters, many of them degenerated.

Concerning the position along the line-of-sight we unfortunately only have limited constraints available for certain sources to fix the three-dimensional structure on all scales. For optically thin radiation, there is no information on the line-of-sight structure. Optically thick radiation comes from a photospheric surface where the opacity exceeds unity. Thus, maps from different wavelengths, with different opacities tracing different radii, can give some constraints on the relative positions of sources along the line-of-sight, but in general it is not possible to fix the absolute positions along the line-of-sight. The more observational data is available the better the full three-dimensional structure could be deduced. We have thus only varied the line-of-sight position in a very limited way. More will be possible when fitting individual molecular lines, because multiple lines from many isotopologues of the same species provide many more surfaces with the opacity exceeding unity, allowing a tomography of the source.

Exploring the complete parameter space of all remaining free parameters is prohibitive in terms of computing time. Thus, we used a hybrid approach to constrain the other parameters mentioned. First, we derive a good guess by varying parameters by hand, and judging the quality visually (χ²-by-eye). This enabled us to fix the exponent, the presence of stellar heating source/star clusters and the displacement in right ascension and declination. We then employed the model optimizer MAGIX (Möller et al. 2013) to iteratively search for the best solution, i.e. the parameter set with the lowest reduced χ² value. For this run, we use the Genetic Algorithm, leaving the central densities and the radii as free parameters. The contribution from the Hii regions (e.g. number electron density) was kept fixed during the entire fitting procedure.

We fit the 874 GHz SMA maps, 480−1280 and 1440−1900 GHz HIFI spectral scans as well as the large scale maps from ATLASGAL at 850 μm, HiGAL SPIRE at 500, 350 and 250 μm and HiGAL PACS at 70 μm.

In total, we have run approximately 2 × 10⁴ different models (multi-core). The runtime for a single run varies between 0.5 and 6.0 h, depending on the density structure, which, e.g. affects the number of final cells in the grid.

4. Analysis and discussion

Within the scope of this paper, we will focus the analysis on the regions for which a wealth of data is available, namely Sgr B2(N), Sgr B2(M), and partly Sgr B2(S). We included Sgr B2(NE) and Sgr B2(deepS) in the model. However, we did not attempt to fit these components, since the available dataset is scarce.

4.1. Large scale dust continuum

To recover the large scale structure of the envelope, as visible in the ATLASGAL and Hi-GAL intensity maps, we superimpose many density components having profiles with varying exponents (see Sect. 3.4). While fitting the large-scale dust continuum maps, we noted the following problem. For a model where the density profile for each component is spherically symmetric and all components are located in the z = 0 plane (Model A), the 70 μm intensity was always overestimated and the 870 μm intensity was underestimated, especially for the two cores Sgr B2(N) and Sgr B2(M).

By moving Sgr B2(N) along the line of sight behind the plane containing Sgr B2(M) as suggested by e.g. Goldsmith et al. (1990) (Model B), the emission from Sgr B2(N) gets attenuated by the envelope of Sgr B2(M) and an improved fit is possible. However, the general trend of the overestimation of the intensity at 70 μm and underestimation of the intensity at 870 μm remains, especially for Sgr B2(M). This is a sign of the dust column density being underestimated and the luminosity being overestimated.

In a first approach, we have reduced this effect by assuming that the envelopes of both, Sgr B2(N) and Sgr B2(M) are elongated along the line-of-sight, i.e. they basically look like a cigar (Model C). We achieved an improved fit by increasing the radii along the line-of-sight by factors of 1.4 and 2.3 for Sgr B2(N) and Sgr B2(M), respectively. At 70 μm the dust becomes optically thick. So by increasing the dust column density along the line-of-sight, we are able to hide the contribution from the stars. On the other hand, this approach increases the intensity at 870 μm, because it increases the total dust column density, which is proportional to the intensity of the optically thin 870 μm emission.

Another solution would be to decrease the luminosity. The luminosity was calculated assuming a single ionized star in each Hii region. However, besides the Lyman continuum emission from the early-type star, additional UV photons could be emitted from accretion shocks in the stellar neighborhood (Cesaroni et al. 2016) which would lead to assigning an earlier spectral type and thus overestimating the total luminosity. To test this possibility, we have multiplied the luminosity by an arbitrary factor of 0.5 (Model D). While doing so, we of course needed to increase the density to preserve the 870 μm fluxes. In this approach, there was no need for the elongated envelope, except for Sgr B2(M), where the envelope is still slightly elongated. However, the observational and theoretical evidence of this Lyman continuum overluminosity seems to be mainly given for B-type stars (Sánchez-Monge et al. 2013; Smith 2014). Thus it is unclear if O-type stars exhibit the same behavior. Furthermore, the effect of dust inside the Hii regions has not been assessed. This will be part of future work, see Sect. 6.

We have kept the number of components in all these models fixed. All parameters of each component are listed in Table B.3. To facilitate the comparison of these different models, we have produced cuts of the intensity along the reversed right ascension axis at the constant declination of the components Sgr B2(NE), Sgr B2(N), Sgr B2(M) and Sgr B2(S). These cuts are shown in Fig. 3. For comparison, we also show the observed and synthetic azimuthally averaged radial profile for all four models of the same components in Fig. A.4. We present the ATLASGAL 870 μm, HiGAL-SPIRE 500, 350 and 250 μm, and the HiGAL-PACS 70 μm from left to right, i.e. the wavelength decreases from left to right. These different models show that there is a degeneracy between the dust density, luminosity and the relative location of the different dust density centers and star clusters along the line of sight. So in summary Model A overestimates the 70 μm intensity for both, Sgr B2(N) and Sgr B2(M) while underestimating the intensity for both regions at 870 μm. Model B only improves the discrepancy for Sgr B2(N), but still gives a bad fit for Sgr B2(M). These two models are thus clearly unfavorable to proceed the analysis with them. Model D provides an improved fit, however it remains unclear to which degree the luminosity could be adjusted. We thus choose to use Model C for the following analysis, since it provides a reasonably good fit for Sgr B2(N) and Sgr B2(M) without assuming any modification of the luminosity.

The resulting large-scale single-dish maps of Model C are shown in Fig. 4. The first row shows the observed maps, the second row the synthetic maps and the third row is a simple difference map between the observed data and synthetics maps. From left to right we have again the same wavelengths as in Fig. 3: ATLASGAL 870 μm, SPIRE 500 μm, SPIRE 350 μm, SPIRE 250 μm, and PACS 70 μm. In the observed maps, the two hot cores Sgr B2(N) and Sgr B2(M) are clearly visible. The extension to the north-east, i.e. Sgr B2(NE), and the one to the south, i.e. Sgr B2(S), are also distinguishable.

In both the observed and the synthetic maps, the intensity of Sgr B2(M) is stronger than Sgr B2(N) except for 870 μm where the opposite is the case. So the model reproduces the general behaviour, but is currently incapable to reproduce the absolute intensities. Especially in Sgr B2(M), the model underestimates the flux at 870 μm and slightly overestimates the flux at 70 μm. The observed maps and the intensity cuts (see Figs. 3 and 4) furthermore show a clear asymmetry, i.e. the profile of the continuum emission has at all wavelengths a slope that is steeper on the east side of the peak than on the west side. To quantify the difference of these slopes, we have fitted the intensity profiles, I ∝ r^{− p}, of Sgr B2(N), Sgr B2(M) and Sgr B2(S) for the wavelengths 870, 500, 350, and 250 μm. For the east side of the peak, we obtain an averaged p of 1.0 and for the west side of the peak we obtain an averaged p of 0.4. Looking at the slopes for the individual sources, we furthermore find that the steepness of the west side of the peak remains unchanged with declination, whereas the steepness of the east side of the peak changes significantly with declination. For Sgr B2(N) and Sgr B2(M), p is on average 1.2, whereas it is 0.7 for Sgr B2(S). This asymmetry is impossible to account for using a single spherically symmetric component. If the slope of the intensity profile, i.e. the column density distribution, is p, then the slope of the density profile should be p + 1( =η). So based on the fitting of the intensity profile, the east wing should have a density exponent of η = 2.0 and the west wing should have a density exponent of η = 1.4.

In the model, we are able to approximate the asymmetry by using a superposition of two components with different exponents for each clump. The east wing is fitted with a density profile having an exponent of η = 2.5, whereas the west wing is fitted with a density profile having a lower exponent of peta = 1.8. These exponents are in good agreement with the ones derived from the intensity profile fitting mentioned above. This clearly shows that our assumption of spherically symmetric clumps as basis limits the modeling effort.

In addition, it is noticeable that the peak of the intensity profiles is shifted to the west in the 870 μm ATLASGAL maps compared to all Hi-GAL maps by ≈0.2 pc, corresponding to ≈5′′ . This is well within the relative pointing error of the ATLASGAL map (pointing accuracy ~4′′ ) and the Hi-GAL maps (pointing accuracy ~1′′ ). Thus we do not attempt to reproduce this shifting ot the peak with our model.

4.2. Small scale dust and free-free continuum

The interferometric maps covering spatial scales up to ~0.25 pc are shown for Sgr B2(M) in Fig. 5 and for Sgr B2(N) in Fig. 6. We refer to these spatial scales as the small-scale structure. In addition to the intensity maps at different wavelengths, we also show the physical setup of Sgr B2(M) in Fig. 7. We plot three density isocontours colored using the dust temperature. To allow us a view inside the model, we cut it open along the (x, y, 0)-plane. The bubbles visible there are the Hii regions that contain no dust. Their dust surface however gets heated by the UV radiation field from the embedded young stellar object. Dust temperatures as high as 600 K are reached.

For the dust cores we find, similar to Qin et al. (2011, for Sgr B2(N)) and Rolffs et al. (2011c, for the hot molecular core G10.47+0.03), a Plummer exponent of η = 5.0 reproducing the small-scale density structure well. For a comparison, we plot the observed and synthetic azimuthally averaged radial profile of each component in Fig. A.5. The parameters of all dust cores are listed in Table B.2. Note that for some of the dust cores, we had to include an internal heating source by placing B-type stars at the center of the core. These could be e.g. stars that do (not yet) show signs of an Hii region. This is also documented in the above mentioned Table B.2. We also had to shift a few of the many Hii regions in Sgr B2(M) behind the dust components to produce a good fit. Please note that at frequencies higher than 100 GHz only one high resolution dataset is available in this modeling procedure. Thus our resulting density distribution is only one possible solution. Other density distributions, e.g. dust cores with a lower exponent and a smaller radius, might yield similar good fits.

However, when comparing Sgr B2(M) and Sgr B2(N) a clear difference is noticeable. Sgr B2(M) appears to be more fragmented compared to Sgr B2(N), which appears rather monolithic. This has already been noticed by Qin et al. (2011). The model presented here allows us, however, to quantify this difference further. We will do this especially in Sect. 4.4. But first we have a look at the stellar population.

4.3. Heating sources

For all four models, we used the same stellar distribution. We account in total for four stellar clusters: Sgr B2(NE), Sgr B2(N), Sgr B2(M), and Sgr B2(S). In Fig. A.3, we plot for each of these clusters a radial distribution histogram of all high-mass stars placed in the model based on observed Hii regions (see Sect. 3.5). We use these histograms to determine by-eye the star cluster radii of 0.8, 0.4, 0.5, and 0.35 pc for Sgr B2(NE), Sgr B2(N), Sgr B2(M), and Sgr B2(S), respectively. Within these radii we sprinkle stars separately for each cluster (see Sect. 3.3). We then treat Sgr B2 as a cluster itself and continue to sprinkle stars everywhere, except in the previously mentioned star clusters, until the calculated total stellar mass of Sgr B2 is reached. We call this part the envelope of Sgr B2. The stellar density distribution is plotted in Fig. 8 and the mass distribution of the stellar mass produced with this approach (including the high-mass stars placed by hand) is shown in Fig. 9. Instead of plotting a histogram, we use a kernel density estimate. A histogram has several disadvantages. It is unsmooth and it depends on the end points as well as on the widths of the bins. By using a smooth kernel we obtain a smooth distribution that is independent of the choice of the end points of the bins and only depends on the chosen bandwidth. For our kernel density estimate a bandwidth of 0.217 M_⊙ has been calculated in linear space.

The results of applying this star sprinkling algorithm are listed in Table 2. For each cluster we list the radius, the number of initial stars, the total number of stars, the enclosed stellar mass, and the luminosity from the initial and all stars, respectively. Where available, we also list luminosities from Lis & Goldsmith (1990) as a reference. Additionally, this table lists the enclosed gas mass for each of the clusters. It is calculated from the dust density distribution assuming a gas-to-dust ratio of 100 (Hildebrand 1983). We calculate the stellar luminosity to mass ratio as 90 L_⊙M_⊙^-1, 71 L_⊙M_⊙^-1, 1176 L_⊙M_⊙^-1, and 120 L_⊙M_⊙^-1 for Sgr B2(NE). Sgr B2(N), Sgr B2(M), and Sgr B2(S), respectively.

In their appendix, Belloche et al. (2013) also estimate the stellar mass of the entire Sgr B2 region by extrapolating the already distributed stellar mass of all stars embedded in observed ultra-compact Hii regions using different IMFs. We apply the same method. However, our total stellar mass is an order of magnitude higher than their result. This difference is explained by the different initial parameters applied and different conversion tables for spectral type to stellar luminosity used. Belloche et al. (2013) calculate a total stellar mass of 675 M_⊙ for the 41 ultra-compact Hii regions they account for. Their lowest mass star has a mass of 11 M_⊙ using the conversion table from Panagia (1973). Using Kroupa’s IMF, they extrapolate a total stellar mass of ~3900 M_⊙ for the mass range of 0.01−120 M_⊙. In our study, we account for more than twice the number of Hii regions, of which the embedded lowest mass star has a mass of 18 M_⊙ using the conversion table from Vacca et al. (1996).

We infer a large population of low- and intermediate mass stars. These low- and intermediate mass stars contain the major fraction of the stellar mass, but the high-mass stars contribute the major fraction of the luminosity. To evaluate the influence of the additional low- and intermediate mass stars, we have performed a simulation excluding these stellar populations. The intensity levels in the large-scale maps are unaffected. We then investigated the dust temperature profiles. Figure 10 shows the density and dust temperature profile along the line-of-sight towards the two hot cores. These line-of-sight profiles are obtained within a 0.5 arcsec beam. Differences in dust temperature on the order of 5−40 K (<10%) are only visible on very small scales towards the cores. The envelope temperature is only affected marginally.

4.4. Mass distribution and star formation efficiency

From the model we obtain the three-dimensional distribution of the stellar and gas mass. We can thus investigate the star formation efficiency (SFE), i.e. the fraction of gas that has been processed into stars. We calculate the SFE using the standard definition (e.g. Myers et al. 1986; Federrath & Klessen 2013): (9)where M_∗ is the stellar mass and M_gas is the gas mass. In Fig. 11 we show the radial distribution of the stellar and gas mass around all four star clusters: Sgr B2(NE), Sgr B2(N), Sgr B2(M), and Sgr B2(S) as well as the derived radially resolved star formation efficiency. For Sgr B2(M) we derive star formation efficiencies of on average 50% , for Sgr B2(N) the derived value is 5% and thus significantly lower. Lada & Lada (2003) suggest that the SFE of a cluster increases with time and can reach a maximum value of typically 30% by the time the cluster emerges from its parental cloud core. This indicates, that Sgr B2(M) has already lost a significant portion of its gas mass through conversion into stars and dispersion. This stellar cluster is emerging from the parental cloud, whereas Sgr B2(N) is still forming stars. This agrees well with the findings by Vogel et al. (1987). For Sgr B2(NE) as well as Sgr B2(S), we calculate rather high average star formation efficiencies of ~30% .

These results have to be taken with caution. The SFE depends on the distributed gas and stellar mass. From the cross cuts in Fig. 3, we think the derived gas masses on large scales are reasonable for all clusters. But while our star sprinkling algorithm will always distribute the same amount of stars if the number of initial stars and their mass is unchanged, the volume within which these stars are placed is set by the star cluster radius. A larger cluster radius effectively decreases the stellar mass density thus causing a lower SFE.

However, to interpret these results, we analyze the relationship between the gas mass of the clusters and their total luminosity. This allows us to study the evolution of the young stellar objects during the phase when they are still embedded in their parental cloud. We follow the same steps outlined in Molinari et al. (2008) and use their evolutionary tracks. In these evolutionary models, the fast accretion phase is stopped when the tracks approach the line where Hii regions are found (at a few thousands solar masses). To match the order of magnitude provided by the massive clusters studied here, we had to arbitrarily extend this line. Please note that these evolutionary models assume that a single massive star is forming, which is clearly not the case for at least Sgr B2(N) and Sgr B2(M). We present all tracks in Fig. 12. The five least massive tracks are the ones presented in Molinari et al. (2008). To obtain the two most massive ones, we had to stop the accretion when the central star reaches 80 and 220 M_⊙, respectively. Based on the instantaneous accretion rates along the tracks, we estimate that in reality the forming star reaches the ZAMS earlier than reaching the Hii line. In particular around 50 M_⊙ for the first track and 75 M_⊙ for the second. The track luminosity where this occurs is around 4.8 × 10⁵ L_⊙ for the first track and around 10⁶ L_⊙ for the second track. If we assume a single massive star is forming, then these tracks suggest that Sgr B2(S) should be compatible with Hii regions (at least one Hii region is observed towards Sgr B2(S)), and Sgr B2(NE) as well as Sgr B2(N) are reaching the ZAMS now. For Sgr B2(M) we do not obtain any useful interpretation. Apart from the fact that a realistic calculation would use a stellar cluster rather than a single star, to reach this relation between gas mass and luminosity in full accretion would require producing a star of a very large mass; or else a star with much larger initial clump mass which would reach the location of M during the ZAMS evolution.

Kruijssen et al. (2015) presented an orbital structure of the Central Molecular Zone (CMZ), i.e. the central few 100 pc of the Milky Way and derived an evolutionary timeline. According to their model, the orbital time between G0.253 (the Brick) and Sgr B2 is Δt = 0.43 Myr. Given the uncertainties in the ages and the age difference, we can only say that the approximate age is compatible with star formation beginning when Sgr B2 was at the position where G0.253 is now. However, the scenario is not able to shed light on the different development stages of the various sources in the region, because on the scales considered by the gaseous streams model all Sgr B2 sources are cospatial.

4.5. Column density map and probability density function

We convert the dust density to gas density by assuming a gas-to-dust ratio of 100 (Hildebrand 1983), and that 73% of the gas is in H₂ (Cox 2000). An H₂ column density map of the whole cloud complex is then obtained by simply summing up the H₂ density distribution along the line-of-sight. The map of the full model with a resolution of 0.1′′ and zoom-ins to Sgr B2(N) and Sgr B2(M) are shown in Fig. A.2. The Hii regions are clearly visible in the zoom-in maps due to their lack of dust.

The propability density function (PDF) of the H₂ column density within molecular clouds is commonly used as a tool to investigate the influence of various competing star formation processes within them (see e.g. Kainulainen et al. 2009; Schneider et al. 2013; Federrath & Klessen 2013).

Unfortunately observational data covering the intermediate scales of Sgr B2 (see Fig. A.6) are lacking. This affects the scales from ~10−20′′ , corresponding to ~0.4−0.8 pc. So while the total flux on these scales is recovered in the Hi-GAL and ATLASGAL maps, its exact distribution is uncertain. We thus refrain from showing the PDF obtained from the H₂ column density map of the model. More observational data is needed.

4.6. Fitting the spectral energy distribution (SED) towards N and M

The fit of the spectral energy distribution towards Sgr B2(N) and Sgr B2(M) obtained with HIFI is shown in Fig. 13. The fluxes from the simulated maps are convolved to the frequency-dependent beam of the 3.5 m Herschel telescope. This kind of fitting thus does not suffer from the necessary convolution to the worst spatial resolution required for modified blackbody fitting. For comparison with previous studies towards Sgr B2(M) and Sgr B2(N) (e.g. Goldsmith et al. 1992; Etxaluze et al. 2013), we have convolved our continuum maps with a beam of 30′′ and extracted the flux towards both sources (see Fig. 14). The agreement between observations and simulation in the ranges, where observational data is available is very good. We note however that the model is not able to reproduce fluxes at wavelengths shorter than 70 μm.

Full three-dimensional continuum radiative transfer modeling has so far been performed for disks, and ρ Oph D (Steinacker et al. 2005). However, this has to our best knowledge not been attempted for Sgr B2, which limits our possibilities of comparison. So in order to compare our results to the work from other authors working on Sgr B2 (e.g. Lis & Goldsmith 1989; Goldsmith et al. 1990; Etxaluze et al. 2013), we have applied the modified black body fitting technique to our synthetic maps.

The intensity I of a blackbody at frequency ν is given by the Planck function (10)where h is the Planck constant, k is the Boltzmann constant, c is the speed of light in vacuum and T is the temperature of the blackbody.

This equation is multiplied by the factor 1−e^{− τ}, where the optical depth τ is given by (11)μ_H2 is the mean molecular weight of hydrogen, m_H is the mass of hydrogen, N(H₂) is the hydrogen column density. The dust opacity κ_ν is determined as a continuous function of frequency ν, by fitting a power-law to the dust opacities given by Ossenkopf & Henning (1994): (12)where κ₀ is the reference dust opacity and χ_d is the gas-to-dust mass ratio.

For this, we have extracted the flux within a beam of 40′′ towards the SPIRE FTS pointing positions, given in Etxaluze et al. (2013), namely α_J2000 = 17^h47^m20.00^s, δ_J2000 = –28°22′17.44″ (Sgr B2(N)) and α_J2000 = 17^h47^m20.30^s, δ_J2000 = –28°23′04.1″ (Sgr B2(M)).

We apply MAGIX (Möller et al. 2013) using first the Genetic Algorithm to find the global best solution, followed by a Levenberg-Marquardt fit to obtain the local best fits for both cores, Sgr B2(M) and Sgr B2(N). The dust temperature T_d, the dust spectral index β and the hydrogen column density N(H₂) where left as free parameters. The mean molecular weight μ_H2 is assumed to be 2.8 (Kauffmann et al. 2008), the gas-to-dust mass ratio χ_d is assumed to be 100 (Hildebrand 1983), the reference dust opacity κ_{500 μm} is 1.773 cm² g^-1 assuming no grain mantles and no coagulation (Ossenkopf & Henning 1994).

We perform a fit covering the frequency range from 450 GHz to 2.87 THz, i.e. excluding the optically thick regime. The resulting SEDs are shown in Fig. 15. Our results are summarized in Table 3. For Sgr B2(N), we obtain a dust temperature T_d of ~31 K, a dust spectral index β of 1.8 and a column density N(H₂) of 2.9 × 10²⁴ cm^-2. For Sgr B2(M) we obtain a dust temperature T_d of ~36 K, a dust spectral index β of 1.9 and a column density N(H₂) of 2.5 × 10²⁴ cm^-2.

These values are the averaged values in a 40′′ beam and are thus the results for the envelopes of Sgr B2(N) and Sgr B2(M). The fitted dust temperature values correspond to radii of 2.0 and 1.3 pc for Sgr B2(N) and Sgr B2(M), respectively (see Fig. 10).

Recently Etxaluze et al. (2013) convolved the Herschel/SPIRE FTS spectral scan maps covering the wavelength range from 194 μm to 671 μm to 40′′ and performed modified black body fits towards Sgr B2(N) and Sgr B2(M). They find spectral indices β of 2.2 and 2.3, dust temperatures T_d of 30 K and 37 K, and column densities N(H₂) of 7 × 10²⁴ cm^-2 and 5 × 10²⁴ cm^-2 for Sgr B2(N) and Sgr B2(M), respectively. Assuming a distance of 8.5 kpc, a mean molecular weight μ_H2 of 2.3, and a dust opacity of κ_{250 μm} of 5.17 cm² g^-1 (Li & Draine 2001), they then estimate dust masses of 2500 M_⊙ and 2300 M_⊙ and luminosities of 1.1 × 10⁶ L_⊙ and 5 × 10⁶ L_⊙ for Sgr B2(N) and Sgr B2(M), respectively.

The derived dust temperatures from our study and the ones from Etxaluze et al. (2013) agree very well. But we obtain lower dust spectral indices for which we see two reasons. On the one hand our modeling setup is limited such that only a single dust species is considered. On the other hand different datasets are considered in both studies. We include the ATLASGAL 870 μm map, as well as the Hi-GAL 500, 350, 250, and 70 μm maps, all cross-calibrated with Planck data and in the case of the Hi-GAL map also cross-calibrated with IRAS data, whereas Etxaluze et al. (2013) only considered Hi-GAL data. This could introduce differences in the fluxes which may result in different dust spectral indices as well as different column densities. So this might also explain why we derive a factor of two difference in the dust column densities.

An important point is the resolution of the observation. We have calculated high resolution (i.e. the pixel size equals 0.1′′ ) synthetic intensity maps covering the same frequency range as the SED plotted in Figs. 13 and 14. We have then convolved these images with decreasing beamsizes from 30′′ down to 0.5′′ . We extract the SEDs towards the Herschel/HIFI positions of Sgr B2(N) and Sgr B2(M) as specified above for each beamsize and also extract the respective τ = 1 surface at each wavelength. This is shown in Fig. 16. From these two SEDs, we can clearly see, that towards the position of Sgr B2(N) the free-free emission becomes optically thicker for smaller beamsizes. This indicates that an Hii region is located along the line-of-sight. Furthermore the dust also becomes optically thick towards Sgr B2(N) with decreasing beamsize, leaving a very narrow frequency range of 100−200 GHz, where the dust is optically thin while still dominating over the free-free emission. The contribution of free-free emission at 100 GHz is in the order of 30%, which is still a significant fraction of the overall emission at this frequency. For Sgr B2(M) the dust remains optically thin for this specific position up to 600 GHz. However, at 100 GHz, the fraction of free-free emission is in the order of 70% for Sgr B2(M), indicating that the free-free emission still dominates at this frequency. Looking at the short-frequency range it is interesting to note that the free-free is partly optically thick for intermediate beamsizes, but turns optically thin for small beamsizes. This shows that there is no Hii region located along the line-of-sight towards the position of the Herschel/HIFI beam on small scales. Given that there are a lot of Hii regions located along the line-of-sight towards Sgr B2(M) this shows that the pointing choice might be crucial for high-resolution observations. It furthermore shows that if instead of a smooth dust density distribution we would assume a clumpy distribution, this clumping would change the photon penetration depth at certain sightlines.

5. Conclusions

We successfully reconstruct a possible three dimensional density distribution of Sgr B2, recovering the continuum structures covering a wide frequency range (ν = 40 GHz−4 THz) on scales from 100 au to 45 pc. We employ the publicly available three-dimensional radiative transfer program RADMC-3D and calculate the dust temperature self-consistently.

• We find that the density field of Sgr B2 can bereasonably well fitted by a superposition of spherical symmetricdensity cores with Plummer-like profiles.

• To reproduce the spectral energy distribution, we locate Sgr B2(N) along the line of sight behind the plane containing Sgr B2(M).

• Sgr B2 comprises a total gas mass of 8.0 × 10⁶ M_⊙ within a diameter of 45 pc. This corresponds to an average gas density of 170 M_⊙ pc^-3. For Sgr B2(N) we find a stellar mass of 2400 M_⊙, a luminosity of 1.8 × 10⁶ L_⊙, a H₂ column density of 2.9.7 × 10²⁴ cm^-2 in a 40′′ beam. For Sgr B2(M) we derive a stellar mass of 20 700 M_⊙, a luminosity of 1.2 × 10⁷ L_⊙, a H₂ column density of 2.5 × 10²⁴ cm^-2 in a 40′′ beam. For Sgr B2(S) we find a stellar mass of 1100 M_⊙, a luminosity of 6.6 × 10⁵  L_⊙, a H₂ column density of 2.2 × 10²⁴ cm^-2 in a 40′′ beam.

• The calculated star formation efficiency is very low for Sgr B2(N), ~5% , and much higher for Sgr B2(M), ~50% . This indicates that most of the gas in Sgr B2(M) has already been converted to stars or dispersed.

6. Outlook

The setup of the density structure and temperature field presented here opens the stage of a plethora of different applications and improvements. We list some of them:

• High resolution, envelope covering maps: with the cur-rent state of telescopes, it is now possible to efficiently maplarger portions of the sky at high resolution. As indicated inFig. A.6, maps with high angular resolution ofSgr B2 are missing. However they are necessary to improve thesetup of the envelope. Currently there is for example only onemap available, that resolves dust structures down to ~0.01 pc towards Sgr B2(N) and Sgr B2(M). To constrainthe dust properties properly, at least another map covering sim-ilar scales at a different wavelength is needed. Studying spatialvariations of the dust properties on scales smaller than 10 ′′ in the entire envelope is impossible with the datasets currently available. However, this could be partly achievedby using the Atacama Compact Array (ACA).

• Physical setup: the physical setup of the model presented here sets the stage for future improvements. By including the possibility to spatially vary the composition of dust, the influence of different compositions can be tested. The treatment of the Hii regions can be improved, eventually leading the way to also model radio recombination lines and thus getting access to the ionized gas content. Furthermore the effect of a clumpy density structure can be tested as well as the influence of dusty Hii regions.

• Molecular line studies: ultimately, one wants to model molecular lines. The setup of the density and temperature structure presented here builds the foundation for such follow-up studies. By keeping the density structure fixed, the remaining main free parameters are the molecular abundance and the velocity field. By modeling a variety of different species, including their isotopologues, different surfaces where the opacity exceeds unity can be sampled, allowing a tomography of the source. Having one model for the whole region, also requires to find a solution of the velocity field that fits the plethora of available molecular line data.

• (M)HD simulations: while the work presented here focuses only on the modeling of a specific source − Sgr B2 − it builds a bridge between theoretical work, e.g. from (M)HD simulations, and observational work.

• Chemical models: combining the radiative transfer modeling efforts with sophisticated chemical modeling, would allow us to fix the molecular abundances.

• Planning observations: the model enables us to examine the full parameter space, finding combinations of e.g. transitions, molecules, wavelengths that would allow us to constrain a degenerate parameter observationally. Thus it is a useful tool to plan new observations. In Fig. 17 for example we show a prediction of high-resolution continuum maps at different wavelengths. They could be observed e.g. with ALMA. These maps indicate that in order to search for fragmentation and determine its level, it is necessary to obtain multi-wavelength observations as there is no single wavelength unaffected by optical thickness effects or free-free continuum contribution.

---

Acknowledgments

We thank the anonymous referee for insightful comments that greatly improved this paper. We furthermore thank C. de Pree for providing the VLA continuum maps presented in Gaume et al. (1995), De Pree et al. (1998). This research is carried out within the Collaborative Research Centre 956, sub-project A6, funded by the Deutsche Forschungsgemeinschaft (DFG). DCL acknowledges support for this work provided by NASA through an award issued by JPL/Caltech. S.L.Q. is partly supported by the NSFC under grant Nos. 11373026, 11433004, by the Top Talents Program of Yunnan Province. This research has made use of NASA’s Astrophysics Data System, Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration et al. 2013), APLpy, an open-source plotting package for Python hosted at http://aplpy.github.com, and the SIMBAD database, operated at CDS, Strasbourg, France. Herschel is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA. HIFI has been designed and built by a consortium of institutes and university departments from across Europe, Canada and the United States under the leadership of SRON Netherlands Institute for Space Research, Groningen, The Netherlands and with major contributions from Germany, France and the US. Consortium members are: Canada: CSA, U.Waterloo; France: CESR, LAB, LERMA, IRAM; Germany: KOSMA, MPIfR, MPS; Ireland, NUI Maynooth; Italy: ASI, IFSI-INAF, Osservatorio Astrofisico di Arcetri-INAF; Netherlands: SRON, TUD; Poland: CAMK, CBK; Spain: Observatorio Astronómico Nacional (IGN), Centro de Astrobiología (CSIC-INTA). Sweden: Chalmers University of Technology − MC2, RSS & GARD; Onsala Space Observatory; Swedish National Space Board, Stockholm University − Stockholm Observatory; Switzerland: ETH Zurich, FHNW; USA: Caltech, JPL, NHSC. The ATLASGAL project is a collaboration between the Max-Planck-Gesellschaft, the European Southern Observatory (ESO) and the Universidad de Chile. It includes projects E-181.C-0885, E-078.F-9040(A), M-079.C-9501(A), M-081.C-9501(A) plus Chilean data. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. The Submillimeter Array is a joint project between the Smithsonian Astrophysical Observatory and the Academia Sinica Institute of Astronomy and Astrophysics and is funded by the Smithsonian Institution and the Academia Sinica.

References

Appendix A: Figures

Appendix B: Tables

Appendix C: Derivation of the physical parameters of Hii regions

We consider Hii regions as Strömgren spheres. They are fully ionized and contain no dust. Here we derive the emission measure, the electron density, and the total flux of ionizing photons for such a region. We assume an electron density distribution that equals n_e within the radius r₀ of the Hii region and is zero elsewhere. We can then calculate the emission measure EM by integrating the number of electrons and number of ions along the line-of-sight. The factor accounts for both sides of the sphere. (C.1)If we assume that the Hii regions only contain hydrogen, it follows that the electron density n_e equals the ion density n_i. Since the electron density n_e is uniform within the Hii region, we can rewrite Eq. (C.1) taking the geometry given in Fig. C.1 into account. where D is the distance to the source, θ₀ = r₀/D is the radius of the Hii region in angular units, and is the angular position in the x−y plane (z = 0). If we want to calculate the maximum emission measure, we set a = 0 and the term vanishes. Introducing the angular source diameter θ_source = 2θ₀, we obtain (C.5)However, neither the emission measure nor the electron density are known. So we need another equation to solve Eq. (C.5). Integrating the specific intensity I_ν over the angular size of the source Ω_source yields the flux density F_ν(C.6)The general solution of the radiative transfer equation for the specific continuum intensity is e.g. given by Rybicki & Lightman (1986) as (C.7)where S_ν is the source function of the Hii region and I₀ is the background intensity. We assume an Hii region with uniform temperature and negligible background temperature, i.e. I_ν(0) ≪ S_ν. The intensity I_ν is then given by (C.8)We have to distinguish two cases for the optical depth τ_ν: (C.9)The source function S_ν for the free-free radiation can be defined using the Planck function B_ν at an electron temperature T_e(C.10)where k is the Boltzman constant, c is the speed of light in vacuum and h is the Planck constant. The second step holds if hν ≪ kT (Rayleigh-Jeans approximation). The optical path length for free-free emission was derived by Oster (1961), an approximation is given by Altenhoff et al. (1960). (C.11)This approximation deviates in the region of interest, i.e. 5 × 10³ ≤ T_e ≤ 1.2 × 10⁴ K and 100 MHz ≤ ν ≤ 35 GHz, less than 10% (Mezger & Henderson 1967).

Putting everything together, and assuming optical thin emission, we obtain (C.12)Usually one fits a Gaussian to the data and extracts the flux density and the full width at half maximum from this fit. However, we determine the total flux in a given circular aperture. We thus calculate the solid angle for the given density distribution as follows. (C.13)Now we obtain the following expression for the emission measure for a circular aperture (C.14)From that we can derive an expression to calculate the electron density (C.15)To derive the flux of ionizing photons, Ṅ_i, we balance the number of recombinations and photoionizations within the Hii region. In a Strömgren sphere, this yields (C.16)where β and β₁ are the rate coefficients for recombinations to all levels and to the ground state, respectively. Thus (β−β₁) provides the recombination coefficient to level 2 or higher. Rubin (1968) approximate the recombination coefficient given by Seaton (1959) for electron temperatures T_e generally found in Hii regions as (C.17)

Including Eq. (C.17) in Eq. (C.16) and solving the integral for a spherical symmetric clump yields: (C.18)Converting the radius r₀ to the angular diameter θ_source of the Hii region and using Eq. (7) yields an expression for the total flux of ionizing photons, given in practical units as (C.19)

All Tables

All Figures