Insights from Synthetic Star-forming Regions. I. Reliable Mock Observations from SPH Simulations

From the simulation output, we know the position of an SPH particle i, its peak density ${\rho }_{i,\mathrm{SPH}}$, its hydrodynamical gas temperature ${T}_{i,\mathrm{SPH}}$, its constant particle mass of ${m}_{i,\mathrm{SPH}}\ ={10}^{-2}\,{M}_{\odot }$ ⁸ and its smoothing length h_i. An SPH particle i has a 3D Gaussian-like profile. The profile strength at position j with distance r_ij from the SPH particle i is defined by the kernel function (for a detailed review, see Springel 2010b):

$$\begin{eqnarray}&&{W}_{{ij}}({r}_{{ij}},{h}_{i})=\displaystyle \frac{1}{{h}_{i}^{3}\pi }\\ \times \,\left\{\begin{array}{lc}1-1.5{({r}_{{ij}}/{h}_{i})}^{2}+0.75{({r}_{{ij}}/{h}_{i})}^{3} & \forall {r}_{{ij}}\lt {h}_{i}\\ \quad 0.25{(2-({r}_{{ij}}/{h}_{i}))}^{3} & \forall \ {h}_{i}\leqslant {r}_{{ij}}\lt 2{h}_{i}.\\ \qquad 0 & \forall {r}_{{ij}}\gt 2{h}_{i}\end{array}\right.\end{eqnarray} \tag{ 1 }$$

The stellar particles in SPH simulations are represented by sink particles (Bate et al. 1995). They form from SPH particles that are bound by the gravitational potential of a forming stellar object and follow the following conditions:

• a.  
  Density threshold. The density threshold is evaluated from the minimal resolvable Jeans mass ${M}_{\mathrm{Jeans}}$ and the isothermal temperature in the simulations. SPH particles, which have higher densities than the threshold density, are tracked and the nearest neighbors are tested by the following conditions.
• b.  
  Jeans limit. The average density and temperature of the group of SPH particles. The resulting mean total mass of the group should be higher than the Jeans mass ${M}_{\mathrm{Jeans}}$.
• c.  
  Potential. The group should be gravitationally bound.
• d.  
  Velocity field. The group should have a negative divergence of the velocity field at the location of the densest particle of the group, which indicates a collapse.
• e.  
  Energy balance. The group of collapsing SPH particles should also have a higher gravitational energy than its rotational kinetic energy to ensure that the collapse cannot be stabilized by angular momentum conservation.

In the D14 SPH simulations, a sink particle forms once more than 50 SPH particles satisfy the above conditions at a given position. In Table 1, we give the number of sink particles at different time-steps in one of the D14 SPH simulations.

Table 1.  Summary of Properties of the run I Simulated Star-forming Region

Time-step	Time	$r({M}_{1/2})$	Sink	${M}_{\mathrm{gas}}$	Feedback
(ID)	(Myr)	(pc)	Particles	$({M}_{\odot })$	on?
024	3.576	8.45	3	9873	no
025	3.725	8.52	4	9825	no
026	3.874	8.60	12	9750	no
027	4.023	8.70	15	9664	no
028	4.172	8.77	17	9593	no
029	4.321	8.83	21	9536	no
030	4.470	8.89	30	9476	no
031	4.619	8.94	32	9414	no
032	4.768	9.01	44	9341	no
042	4.917	9.08	48	8833	yes
052	5.066	9.27	48	8643	yes
062	5.215	9.58	53	8598	yes
072	5.364	9.91	61	8495	yes
082	5.513	10.30	73	8377	yes
092	5.662	10.79	85	8187	yes
102	5.811	11.30	91	7987	yes
112	5.960	11.84	98	7722	yes
122	6.109	12.37	100	7500	yes
132	6.258	12.94	108	7254	yes
142	6.407	13.50	115	6996	yes
152	6.556	14.09	122	6687	yes
162	6.705	14.69	125	6409	yes
172	6.854	15.33	128	6105	yes

Download table as:  ASCIITypeset image

From the minimum resolved Jeans mass ${M}_{\mathrm{Jeans}}$, a respective length scale, called the accretion radius ${r}_{\mathrm{acc}}$, is calculated, which is comparable to the Jeans length ${\lambda }_{\mathrm{Jeans}}$. SPH particles are accreted by a sink particle when they are within the accretion radius ${r}_{\mathrm{acc}}$, are gravitationally bound by the sink particle's potential and have low enough rotational energies. The sink particles⁹ grow over time through accretion of many SPH particles.

High-mass stellar feedback is switched on in the D14 SPH simulations as soon as three high-mass sink particles above 20 ${M}_{\odot }$ have formed.¹⁰ The simulation ends before the first supernova of a high-mass sink particle would presumably go off.

In this work, we explore run I of the D14 SPH simulations an intermediate mass cloud with a total gas mass of 10⁴ ${M}_{\odot }$ and initially 10⁶ SPH particles. We choose this run because it is of relatively high resolution with a sink accretion radius of ${r}_{\mathrm{acc}}=0.005\,\,\mathrm{pc}$. Individual stellar objects in the form of sink particles with masses between 0.5 ${M}_{\odot }$ and 70 ${M}_{\odot }$ are produced and can be resolved, and therefore can be interpreted as single stars. The first stars formed after about 3.6 $\mathrm{Myr}$ after the start of the simulation (see also Section 2.1.6). The stars in the simulation grow through accretion (see also Section 2.1.7) and ensure that the stellar population aquires its own initial mass function (IMF) without imposing an initial IMF to the simulation (see Dale et al. 2012).

We choose run I because we aim to produce synthetic observations that are as close to real star-forming regions as possible. Therefore, for our studies, we limit ourselves to the coupled feedback scenario: ionization and winds. The ionizing sources in the simulations emit a flux of ionized photons ranging from 1.3 × 10⁴⁸ s⁻¹ to 1.3 × 10⁴⁹ s⁻¹. They also emit winds¹¹ with a speed ranging from 1780 Km s⁻¹ to 3125 Km s⁻¹ and a respective mass-loss rate between 3.1 × 10⁻⁷ ${M}_{\odot }$ yr⁻¹ and 2.8 ×10⁻⁶ ${M}_{\odot }$ yr⁻¹. For a detailed description of the high-mass feedback formalism, see Dale (2015), Dale et al. (2007, 2013b, 2012), and Dale & Bonnell (2008).

In run I, the first stars formed after about 3.6 $\mathrm{Myr}$, high-mass stellar feedback is switched on at about 4.8 $\mathrm{Myr}$ and the simulation is stopped at about 7 $\mathrm{Myr}$ after the start of the simulation. The simulation is stopped before the first supernova¹² would presumably go off. We select 23 time-steps with a constant step size Δt = 149,000 years, starting once the first star has formed. At every time-step ${t}_{\mathrm{step}}$ and for every sink particle, we keep track of the physical properties, such as the sink particle mass ${M}_{\mathrm{sink}}({t}_{\mathrm{step}})$, the age of the sink particle ${t}_{\mathrm{age}}({t}_{\mathrm{step}})$, and the accretion rate ${\dot{M}_{\mathrm{acc}}}({t}_{\mathrm{step}})$. For more details see Table 1.

The sink particle mass is defined as the sum over all ${N}_{\mathrm{acc}}$ accreted SPH particles i:

$$\begin{eqnarray}&&{M}_{\mathrm{sink}}({t}_{\mathrm{step}})=\underset{i=1}{\overset{{N}_{\mathrm{acc}}}{\displaystyle \sum }}{m}_{i,\mathrm{SPH}}.\end{eqnarray} \tag{ 2 }$$

Hereafter, we define the age of a sink particle as

$$\begin{eqnarray}&&{t}_{\mathrm{age}}({t}_{\mathrm{step}})={t}_{\mathrm{step}}-{t}_{\mathrm{appear}}+{\rm{\Delta }}t/2,\end{eqnarray} \tag{ 3 }$$

where ${t}_{\mathrm{appear}}$ is the time of the time-step when a sink particle first forms (binds more then 50 SPH particles). We add the time ${\rm{\Delta }}t/2$ to the age because the sink particle may have formed at any time between the time-step when it appeared and the previous one, and therefore, on average, we assume that it formed half way in between the two.

We consider a sink particle to be accreting between two different time-steps, when

$$\begin{eqnarray}&&{\rm{accretion}}{| }_{{t}_{\mathrm{step}}\to {t}_{\mathrm{step}+1}}={M}_{\mathrm{sink}}({t}_{\mathrm{step}+1})\gt {M}_{\mathrm{sink}}({t}_{\mathrm{step}}).\end{eqnarray} \tag{ 4 }$$

The accretion definition becomes important when mapping the SPH simulation to the Voronoi mesh as described in Section 3. We define the accretion rate at a time ${t}_{\mathrm{step}}$ as the mass-gain rate between two time-steps:

$$\begin{eqnarray}&&{\dot{M}_{\mathrm{acc}}}({t}_{\mathrm{step}})=\displaystyle \frac{{M}_{\mathrm{sink}}({t}_{\mathrm{step}+1})-{M}_{\mathrm{sink}}({t}_{\mathrm{step}-1})}{2{\rm{\Delta }}t}.\end{eqnarray} \tag{ 5 }$$

While the D14 SPH simulations have high-mass stellar feedback mechanisms implemented, the simulations do not take into account the dust heating from the stars. In D14 SPH simulations, the external heating through the diffuse Galactic radiation field is approximated by the Larson (2005) EOS¹³ rather than just an isothermal approximation. The EOS contains the following regimes (following Equation (1) and Figure 1 of Dale et al. 2012):

• a.  
  line cooling processes

  $$\begin{eqnarray}&&P\sim {\rho }^{3/4}\\ &&\forall \ \rho \lt 5.5\times {10}^{-19}\,{\rm{g}}\,{\mathrm{cm}}^{-3},\end{eqnarray} \tag{ 6 }$$
• b.  
  dust cooling processes

  $$\begin{eqnarray}&&P\sim {\rho }^{1.0}\\ &&\forall 5.5\times {10}^{-19}\,{\rm{g}}\,{\mathrm{cm}}^{-3}\lt \rho \lt 5.5\times {10}^{-15}\,{\rm{g}}\,{\mathrm{cm}}^{-3},\end{eqnarray} \tag{ 7 }$$

where P is the gas pressure and ρ is the gas density. Since in the ideal gas approximation (see Feynman 1977), the pressure P always follows

$$\begin{eqnarray}&&P={{Nk}}_{{\rm{B}}}T/V=\rho {k}_{{\rm{B}}}T/{m}_{p}\sim \rho T,\end{eqnarray} \tag{ 8 }$$

with gas particle number N, the Boltzmann constant k_B, the gas temperature T, the gas particle mass m_p, and the gas Volume V, we can also express the line cooling term from Equation (6) as a power law of density ρ and temperature T:

$$\begin{eqnarray}&&P\sim \rho T\sim {\rho }^{3/4}\to \rho \sim {T}^{-4}.\end{eqnarray} \tag{ 9 }$$

At high densities, dust cooling is assumed to dominate and the EOS becomes isothermal with ${T}_{\mathrm{iso}}=7.5\,{\rm{K}}$. The densities in large parts of the cloud are low enough that the gas and dust should not be well coupled thermally.

In Figure 2(a), we have plotted temperature versus density for time-step 122 (at 6.109 $\mathrm{Myr}$) of run I 1.341 $\mathrm{Myr}$ after switching on high-mass stellar feedback (see Table 1). Points regardless of color (blue, red, and black points together) represent all 10⁶ SPH particles of run I. The red and black points represent SPH particles that lie within our selected box of 30 pc width. The red points are (partly or fully) ionized particles (see with Dale et al. 2007, 2012) within the box. The horizontal line at 10⁴ K shows all the fully ionized SPH particles. The other red points are partly ionized and lie within the selected box. The black points are neutral SPH particles in the selected box. The state of these neutral particles is described by the EOS, which is describing the dependency of density and temperature (diagonal line: Equations (6), (8) and (9); horizontal: Equation (7)).

Zoom In Zoom Out Reset image size

In this work, we will use a dust radiative transfer code; therefore, we will only use SPH particles in the box that are completely neutral since dust gets destroyed during the ionization process (Diaz-Miller et al. 1998). Hereafter, we call the clipping for size and only neutral SPH particles, resulting in the black particles in Figure 2(a), method c1.

In Figure 2(b), we can see that the SPH temperature increases with radius from the center of the simulation. Hot ionized particles (cyan) move to larger radii, but SPH particles that still follow the EOS also have higher temperatures at larger radii. This is a result of the EOS, since the density decrease with radius in centrally condensed clouds. These can also be seen in Figure 2(b) and the 2D temperature plot of Figure 2(c), where we only plot the black sample of Figures 2(a) and (b). The maximum temperature at this time-step at the boundary of the box lies around 400 K.

The SPH temperature of neutral particles in the box is lower than the typical dust sublimation temperature of 1600 K (e.g., Whitney et al. 2003), but much higher than the expected ambient dust temperature around a star-forming region of about 18 K (see Paper II). We explore further temperature cuts of the line cooling part of the EOS. Hereafter, we call the clipping in size using only neutral particles, and removing particles above a certain threshold temperature clipping, method c2. Note that method c2 is the same as c1, but additionally particles above a certain threshold temperature are removed.

In Figure 3, we show the temperature histogram of neutral particles with a box of 30 pc of the four time-steps 024 (3.576 $\mathrm{Myr}$), 032 (4.768 $\mathrm{Myr}$), 072 (5.364 $\mathrm{Myr}$), and 122 (6.109 $\mathrm{Myr}$). We can see that the temperature peak shifts from $T=60\,{\rm{K}}$ to about $T=35\,{\rm{K}}$ once the high-mass stellar feedback is turned on after time-step 32 (4.768 $\mathrm{Myr}$). This is due to the sweeping up of the low-density neutral gas into dense, cool clumps.

Zoom In Zoom Out Reset image size

Therefore, a clipping, and hence removing of the SPH particles above a certain temperature (e.g., 60 K) in clipping method c2, is not helpful, since we would lose too many particles for certain time-steps. Also, a threshold temperature below the sublimation is hard to pick, since it remains unclear when exactly the temperature of the dust decouples from the temperature of the gas. Hereafter, we will only remove non-neutral SPH particles from the simulation, as in clipping method c1, and keep SPH particles at all temperatures within the box, which is always below the sublimation temperature.

Note that such a steep increase in temperature with radius is not observed in the dust species of real star-forming regions. Hence, implementing the hydrodynamical temperature, which is a result of the line cooling part (Equation (6)) of the EOS, and therefore gas phase of the region, should happen with caution. We will address this in detail in Section 6, where we present a better treatment of the dust temperature.

A first order approach is to put sites at the positions of the SPH particles of the clipped sample c1 (for details about the sample selection, see Section 2.2). Hereafter, the site method s1 is defined by sites at SPH positions alone. While this works rather well for isolated sink particles (see sink particle 110753 in Figure 4(a)), very large high-density Voronoi cells are created close to the ionizing bubble (see 11640 and 197247 in Figure 4(a)). These large high-density cells are an artifact of the SPH to Voronoi mesh mapping and were not part of the initial simulations. Artificial high-density regions create overestimated fluxes in the radiative transfer because now very large high-density regions lie close to very luminous objects.

In Figure 5, we show the radial peak density distribution of SPH particles close to sink particles for an accreting source in a small cluster (11640), an isolated accreting source (110753) and a source in the ionizing bubble (115358), which lost all its material due to a close-by ionizing source; all profiles of all sources of time-step 122 (6.106 Myr) are displayed in the online extension of Figure 5. Red vertical lines highlight the distance of nearby accreting sink particles, while blue lines highlight nearby sink particles in the ionizing bubble. We can see that accreting sources can cause a radial trend in the distribution of the SPH particles, which can be seen very clearly for sink particle 110753 directly in Figure 5 (middle), but also for the density peaks in the density profile located 1 pc from sink particle 11640 in Figure 5 (middle) and of sink particle 110753 in Figure 5 (left).

Zoom In Zoom Out Reset image size

View extended figure (2 panels)

Ionizing sources create (large) voids in the profile (for 115358 about 0.6 pc). Physically, this void represents ionized low-density bubbles around the ionizing stars. Since there is a low-density region around the ionizing sink particles, very few SPH particles are located close to the other sink particles in the ionizing bubble. For example, the closest SPH particle to sink particle 115358 is about 0.6 pc away.

The Voronoi cells around the sink particles are created by distant SPH particles beyond the bubble that have a higher density than the bubble. When ionizing sources are present in the particle-based simulation, artificial large density cells can be avoided by putting additional sites at the positions of the sink particles. Hereafter, we will refer to this method as site method s2, which combines sites at the position of the SPH particles and the sink particles.

When adding sites at the position of the sink particles, we can reproduce the low-density bubbles around the ionizing sink particles. Hence, the size of high-density cells reaching inside the bubbles is truncated. In Figure 4(b), we show this effect for sink particle 11640 and sink particle 197247, which have ionizing bubbles close-by.

When inspecting Figure 4(b), we see that there are still some large high-density cells relatively close to the accreting sink particles, which is partly due to the lack of resolution of the SPH simulation. Adding additional circumstellar sites around the accreting sink particles increases the resolution around the accreting sink particles and minimizes the artificial high-density cells as can be seen in Figure 4(c). Hereafter, site method s3 combines sites at the position of SPH particles, the sink particle position, and the circumstellar sites around accreting sink particles.

On average, the distance between an accreting sink particle and the closest SPH particle is 0.01 pc, which is about 2000 au as can be seen for the accreting sink particles 11640 and 110753 in Figure 5 (left, middle). Concretely, in the method s3 setup, the position of the sites is constructed using a constant probability distribution function (PDF; for more details, see Appendix B and Press et al. 1992) in log-space of radius

$$\begin{eqnarray}&&r={10}^{\xi [{\mathrm{log}}_{10}{r}_{\max }-{\mathrm{log}}_{10}{r}_{\min }]+{\mathrm{log}}_{10}{r}_{\min }},\end{eqnarray} \tag{ 10 }$$

with ${r}_{\min }={10}^{-4}\,\mathrm{pc}$ and ${r}_{\max }={10}^{-1}\,\mathrm{pc}$. The angles θ and ϕ are spaced isotropically within the sphere:

$$\begin{eqnarray}&&\theta =\arccos (2\xi -1)\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&\phi =2\pi \xi .\end{eqnarray} \tag{ 12 }$$

Once we set the circumstellar sites of method s3 in this way, we can approximate the density gradient of the envelope more smoothly. This will be discussed in more detail in Section 5.1.

In this section, we already explained in detail different methods of how to distribute the sites in a Voronoi mesh with a particle-based code as input. Site distribution method s3, with sites at the neutral SPH positions, at the sink particle positions and circumstellar sites around the accreting sink particles, can increase the resolution and remove artifacts from the particle to mesh transition. This is why we highly recommend distributing sites in this manner. We found that the number of circumstellar sites ${N}_{\mathrm{circ}}\approx 1000$ is a good estimate to increase resolution around the sink particles but not create too many cells, which might slow down computation for the steps described in the next sections. In the remainder of this paper, we explicitly use the site distribution method s3.

To evaluate the density at a new site j at the position of the sink particle or at a circumstellar site around a sink particle, we compute the SPH density function (for a detailed review, see Springel 2010b) for the N closest SPH particles,

$$\begin{eqnarray}&&{\rho }_{j,{\rm{p}}1}=\underset{i=1}{\overset{N}{\displaystyle \sum }}{m}_{i,\mathrm{SPH}}{W}_{{ij}}({r}_{{ij}},{h}_{i}),\end{eqnarray} \tag{ 13 }$$

where ${m}_{i,\mathrm{SPH}}$ is the constant SPH particle mass with the value ${m}_{i,\mathrm{SPH}}={10}^{-2}\,{M}_{\odot }$ in run I of the D14 SPH simulations and the kernel function ${W}_{{ij}}({r}_{{ij}},{h}_{i})$ from Equation (1).

We estimate the temperature at a new site j by weighting (Press et al. 1992) the SPH particle temperature ${T}_{i,\mathrm{SPH}}$ of the closest N particles with distance¹⁵ r_ij

$$\begin{eqnarray}&&{T}_{j,{\rm{p}}1}=\displaystyle \frac{{\displaystyle \sum }_{i=1}^{N}{T}_{i,\mathrm{SPH}}\,{m}_{i,\mathrm{SPH}}{W}_{{ij}}({r}_{{ij}},{h}_{i})}{{\rho }_{j,{\rm{p}}1}}\end{eqnarray} \tag{ 14 }$$

with the help of Equation (13). We explored the effect of setting N to different values and found that setting $N\gt {N}_{\mathrm{neighbor}}$ produces a good approximation¹⁶ of the exact solution of the density. Hereafter, we call this property evaluation method p1.

In Figure 6(a), we plot the density evaluation with property method p1 and site method s3 (note, the different scale than in Figure 4(c)). We can see that even with the site distribution method s3, we still have artifacts of large high-density cells as already mentioned in Section 3.2. These large high-density cells are due to the density overestimation when using the peak density ${\rho }_{i,\mathrm{SPH}}$ at the SPH particles, which lie in large cells.

Zoom In Zoom Out Reset image size

When mapping the density with method p1, we recover a total mass of 139945 ${M}_{\odot }$ which is several orders of magnitude higher than the actual mass of about 7500 ${M}_{\odot }$. Therefore, we explore other techniques.

In the community of SPH modelers, the density is sometimes approximated in Voronoi cells using an average density $\bar{\rho }={m}_{\mathrm{SPH}}/{V}_{\mathrm{cell}}$ (D. Hubber, private communication). We do not use this technique, since it can only be used to evaluate the density in cells containing SPH particles. When adding new sites, this approximation becomes untrue because the volume ${V}_{\mathrm{cell}}$ shrinks, while ${m}_{\mathrm{SPH}}$ remains constant by definition. Also, the density evaluation of cells, which contain the new sites is not covered by that technique.

Thus, we virtually split the mass of an SPH particle i according to the 3D Gaussian-shaped SPH kernel function (see Equation (1)). We repeat this for all SPH particles. Through simple counting, we reevaluate the mass in the cells of the sites j of the existing SPH particle sites plus the newly added sites of method s3. Hereafter, we call this property evaluation method p2.

To space the split SPH particle k with respect to the SPH particle i again, we use a linear spaced PDF (as described in detail in Appendix B) and solve for the positions of the split particles r_ik. We distribute the split SPH particles isotropically (see Equations (11) and (12)). We convert all the ${N}_{\mathrm{split}}^{{\rm{p}}2}$ split SPH particles to Cartesian coordinates and scale the distance by the respective smoothing length h_i of the SPH particle which got split, since we set ${h}_{i}=1$ for the calculation above. The split SPH particles do not necessarily fall in the cell of the SPH particles they split from. By searching for the nearest neighbor sites (Press et al. 1992) to each split SPH particle, we can now reassess the mass of the Voronoi cells, the split SPH particles fall into. We then calculate the new density ${\rho }_{j,{\rm{p}}2}$ at a cell site j as

$$\begin{eqnarray}&&{\rho }_{j,{\rm{p}}2}=\underset{k=1}{\overset{{n}_{j}}{\displaystyle \sum }}{m}_{\mathrm{SPH},{\rm{k}}}/{V}_{j},\end{eqnarray} \tag{ 15 }$$

with the split SPH particle mass ${m}_{\mathrm{SPH},{\rm{k}}}={m}_{\mathrm{SPH},{\rm{i}}}/{N}_{\mathrm{split}}^{{\rm{p}}2}$ and the number n_j set to the number of split SPH particles found in a cell with the volume V_j. We found that splitting by ${N}_{\mathrm{split}}^{{\rm{p}}2}=100$ is a good trade-off between finding a smooth density distribution and computational efficiency.

Since the probability of accumulating many split SPH particles decreases with decreasing cell size, we introduce a threshold of ${N}_{\mathrm{threshold}}^{{\rm{p}}2}=10$ split SPH particles. For cells that collect fewer than ${N}_{\mathrm{threshold}}^{{\rm{p}}2}$, we evaluate ${\rho }_{j,{\rm{p}}1}$ from Equation (13).

In Figure 6(b), we plot the density evaluated with property method p2. When comparing to Figure 6(a), one can see that the transition between cells is smoother and artificial high-density Voronoi cells are substantially fewer.

To evaluate the temperature, we simply keep track from which original SPH particle of mass ${m}_{\mathrm{SPH},{\rm{i}}}$ a split SPH particles of mass ${m}_{\mathrm{SPH},{\rm{k}}}$ comes from and what temperature T_i it had, we can weigh the temperature accordingly by the density:

$$\begin{eqnarray}&&{T}_{j,{\rm{p}}2}=\underset{k=1}{\overset{{n}_{j}}{\displaystyle \sum }}{T}_{i}({m}_{\mathrm{SPH},{\rm{k}}})/{n}_{j}.\end{eqnarray} \tag{ 16 }$$

For the cells with ${n}_{j}\lt {N}_{\mathrm{threshold}}^{{\rm{p}}2}$, we calculate the temperatures as in method p1.

When mapping the density with method p2, we recover a total mass of 7530 ${M}_{\odot }$ relatively close to the real mass of 7500 ${M}_{\odot }$.

Another possibility is to evaluate the density by sampling ${N}_{\mathrm{random}}$ random positions l in the grid and calculating the density from the kernel distribution (as in method p1) of the N closest SPH particles to the positions l.

By averaging¹⁷ the density of the random points l at n_j positions in cell j, we get the total density ${\rho }_{j,{\rm{p}}3}$ of method p3

$$\begin{eqnarray}&&{\rho }_{j,{\rm{p}}3}=\underset{l=1}{\overset{{n}_{j}}{\displaystyle \sum }}{\rho }_{l,{\rm{p}}1}/{n}_{j}.\end{eqnarray} \tag{ 17 }$$

Typically, we set ${N}_{\mathrm{random}}\approx 2e7$ for ${N}_{\mathrm{SPH},{\rm{i}}}\approx 7e5$ neutral SPH particles in the simulation time-step and ${N}_{\mathrm{split}}^{{\rm{p}}3}=30$ as above. To ensure that every cell j has at least ${N}_{\mathrm{threshold}}^{{\rm{p}}3}=20$, we add points to under-filled cells until ${n}_{j}={N}_{\mathrm{threshold}}^{{\rm{p}}3}$.

Figure 6(c) shows the property evaluation method p3. Density gradients are smoother and no artificial high-density cells remain.

The temperature of property method p3 is a weighted version of Equation (17) using Equation (13)

$$\begin{eqnarray}&&{T}_{j,{\rm{p}}3}=\displaystyle \frac{1}{{n}_{j}}\underset{l=1}{\overset{{n}_{j}}{\displaystyle \sum }}\displaystyle \frac{{\displaystyle \sum }_{i=1}^{N}{T}_{i,\mathrm{SPH}}{m}_{i,\mathrm{SPH}}{W}_{{il}}({r}_{{il}},{h}_{i})}{{\rho }_{l,{\rm{p}}1}},\end{eqnarray} \tag{ 18 }$$

again $N\gt {N}_{\mathrm{neighbor}}$ as in method p1.

When mapping the density with method p1, we recover a total mass of 7495 ${M}_{\odot }$ extremely close to the actual mass of 7500 ${M}_{\odot }$.

In Figure 6, we show different methods of mapping the SPH density and temperature onto the Voronoi mesh.

We found that kernel evaluation method p1 is efficient and reliable in small cells. The solving of the kernel function using the $N\gt {N}_{\mathrm{neighbor}}$ nearest neighbors is sufficient to approximate the properties at the new position. But method p1 creates artifacts for large cells.

The splitting kernel method p2 is less computationally efficient but produces a smooth density or temperature gradient between large cells. For small cells, this method breaks down, since the probability that a small cell gets "hit" by split particles is small. Therefore, the introduced threshold ${N}_{\mathrm{threshold}}^{{\rm{p}}2}=10$ is important. For such small cells, we use the method p1. We found that virtually splitting the SPH particle by ${N}_{\mathrm{split}}^{{\rm{p}}2}=100$ produces a solid result while being tolerably efficient in time.

Method p3 removes the statistical sampling issues of the small cells (compared to method p2) by adding random sampling points within the cell as long as ${N}_{\mathrm{threshold}}^{{\rm{p}}3}$ is reached. We found that a smooth gradient in temperature or density can be reached when setting the number of points to sample to ${N}_{\mathrm{random}}\approx 2\times {10}^{7}$ and the threshold to ${N}_{\mathrm{threshold}}^{{\rm{p}}3}=20$. The trade-off is that this method is less efficient computationally when comparing with method p2, but it is fast enough for the purposes of our work.

Below, we list the total mass within the simulation box of time-step 122 (6.109 $\mathrm{Myr}$), recovered from the density mapping for every method present in this section:

$$\begin{eqnarray}&&{m}_{{\rm{p}}1}^{\mathrm{tot}}\approx 139945\,\,{M}_{\odot },\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&{m}_{{\rm{p}}2}^{\mathrm{tot}}\approx 7530\,\,{M}_{\odot },\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&{m}_{{\rm{p}}3}^{\mathrm{tot}}\approx 7495\,\,{M}_{\odot }.\end{eqnarray} \tag{ 21 }$$

The actual mass within the box of the SPH simulation ${m}_{\mathrm{SPH}}^{\mathrm{tot}}=7500\,\,{M}_{\odot }$ (see Table 1). In respect of mass conservation, we note that method p1 is completely unreliable. The small overestimate of method p2 results from the threshold mechanism.¹⁸ In method p2, when there are less than ${N}_{\mathrm{threshold}}^{{\rm{p}}2}$ split SPH particles in the cell, only the density at the position of the site is estimated. This leads to an over-prediction of the cell's density. For this reason, method p1 and method p2 are faulty by definition, since they introduce too high densities to certain cells. Therefore, we suggest method p3 when mapping SPH distributions onto a Voronoi tessellation, since this overestimation is reduced (error below 1%) by sampling more positions in small cells. In Figure 8, we show how well method p3 works in recovering a density gradient from the initial simulation. Blue points represent the simulated peak density of the SPH particles and the black points represent the p3 circumstellar sites. Figure 8 shows that the technique succeeds to recover the density structure over two orders of magnitude over less than 0.1 pc.

Zoom In Zoom Out Reset image size

The Ulrich envelope profile (Ulrich 1976) is a rotational flattened profile and follows the physical process of accretion. The profile ${\rho }_{\mathrm{Ulrich}}(r,\theta )$ is described in spherical coordinates and is dependent on the infall rate of material ${\dot{M}}_{\mathrm{env}}$, the mass of the stellar object M_*, and the centrifugal radius R_c:

$$\begin{eqnarray}{\rho }_{\mathrm{Ulrich}}(r,\theta ) & = & {\rho }_{u}^{\mathrm{env}}{\left(\displaystyle \frac{r}{{R}_{c}}\right)}^{-3/2}{\left(1+\displaystyle \frac{\cos \theta }{\cos {\theta }_{0}}\right)}^{-1/2}\\ & & \times {\left(\displaystyle \frac{\cos \theta }{\cos {\theta }_{0}}+\displaystyle \frac{2{R}_{c}{\cos }^{2}{\theta }_{0}}{r}\right)}^{-1}\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&{\rho }_{u}^{\mathrm{env}}=\displaystyle \frac{{\dot{M}}_{\mathrm{env}}}{4\pi {({{GM}}_{* }{R}_{c}^{3})}^{-1/2}}.\end{eqnarray} \tag{ 24 }$$

The angle ${\theta }_{0}$ is given by the streamline equation:

$$\begin{eqnarray}&&0={\cos }^{3}{\theta }_{0}+\cos {\theta }_{0}\left(\displaystyle \frac{r}{{R}_{c}}-1\right)-\cos \theta \left(\displaystyle \frac{r}{{R}_{c}}\right).\end{eqnarray} \tag{ 25 }$$

Hereafter, we call the Ulrich envelope profile¹⁹ method e1. We calculate the infall rate ${\dot{M}}_{\mathrm{env}}$ from the accretion rate (see Equation (5)). As can be seen in Figure 8, we found that an Ulrich profile (left panel, purple lines) derived from the parameters of sink particle 110753 seamlessly matches the accretion stream inwards but also underestimates it at larger radii. Note that the singularity at the centrifugal radius R_c (when $\cos \theta ={\cos }^{3}{\theta }_{0};$ see Equation (25)) at the mid-plane (for θ = 90 °; see Equation (23)) is caused by the pile-up of material as a result of angular momentum conservation. In reality, the material would be distributed radially to form a circumstellar disk.

The Ulrich description of the circumstellar disk (i.e., the singularity) is not physically realistic. Therefore, for the density profile, we ignore the singularity at the mid-plane and later add a more advanced description of the circumstellar disk in Section 5.4. In this method, hereafter named e2, we suppress²⁰ the Ulrich singularity by adjusting the Ulrich density profile of Equation (23):

$$\begin{eqnarray}{\rho }_{\mathrm{Ulrich}}^{\mathrm{no}\ \mathrm{singularity}}(r,\theta ) & = & {\rho }_{u}^{\mathrm{env}}{\left(\displaystyle \frac{r}{{R}_{c}}\right)}^{-3/2}{\left(1+\displaystyle \frac{\cos \theta }{\cos {\theta }_{0}}\right)}^{-1/2}\\ & & \times {\left(1+2\displaystyle \frac{{R}_{c}}{r}\right)}^{-1}.\end{eqnarray} \tag{ 26 }$$

In Figure 8 (middle), we present the rotationally flattened Ulrich profile with the singularity removed (green). We can see that the Ulrich distribution matches the sites at the inner rim of the envelope resolved by the SPH simulation, but underestimates the resolved SPH envelope profile at larger radii for object 110753.

As an alternative, and in order to understand how important the assumption of the envelope density profile is for the final results, we also consider the case where we use envelopes with spherically symmetric power-law profiles rather than the Ulrich collapse model (referred to as method e3). In the right panel of Figure 8, we show four different radial power-law profiles (diagonal dashed lines) with exponents $\gamma =[-0.5,-1.0,-1.5,-2.0]$ for the density function²¹

$$\begin{eqnarray}&&\rho (r)={\rho }_{0}{(r/{r}_{0})}^{\gamma }.\end{eqnarray} \tag{ 27 }$$

The orange solid diagonal line is the log-scaled least-squares fit (Press et al. 1992) of the density distribution of the nearby SPH particles. The log-scaled least-square fit requires us to fit the following function.

$$\begin{eqnarray}&&\,{\mathrm{log}}_{10}\,\rho (r)=\,{\mathrm{log}}_{10}\,{\rho }_{0}+\gamma \,{\mathrm{log}}_{10}(r/{r}_{0}).\end{eqnarray} \tag{ 28 }$$

For this accreting sink, the best fit of $\gamma =-1.25$ is a typical value for a protostar.

Note that a fitted power-law profile in density typically produces a better fit in terms of absolute scaling compared to the Ulrich profiles. However, spherical symmetric power-law density profiles with only one slope are not entirely physically motivated, whereas the Ulrich models are self-consistent with respect to the infall rate.

The Ulrich description (method e1) is a physically self-consistent envelope-infall model and is therefore consistent with the general physical picture, when improving the resolution on small scales by extrapolating inwards. However, the singularity at the centrifugal radius is unrealistic. Therefore, due to the singularity at the mid-plane, we do not use the standard Ulrich envelope (method e1).

However, we will use the Ulrich description without singularity e2 and add a circumstellar disk later in the radiative transfer setup in Section 5.5. For the Ulrich description, we set the centrifugal radius R_c, hence the outer rim of the circumstellar disk, to ${R}_{c}=500\,\,\mathrm{au}$. The choice of the centrifugal radius is dependent on the chosen outer disk radius. 500 au is a realistic intermediate value of typical values between 100 au and 1000 au. Moreover, the disk radius needs to be smaller than the accretion radius from the simulation (${r}_{\mathrm{acc}}=0.005\,\,\mathrm{pc}\approx 1000\,\,\mathrm{au}$), because otherwise the disk would have been resolved by the simulation. We calculate the density ${\rho }_{u}^{\mathrm{env}}$ (see Equation (23)) for all Ulrich envelope fits. In Figure 9 (top panel), we show e2 rotationally flattened envelope profiles for the three accreting sink particles 11640, 110753, and 197247 with the corresponding values for ${\rho }_{u}^{\mathrm{env}}$ and ${\dot{M}}_{\mathrm{env}}$. To see all profiles for time-step 122 (6.109 $\mathrm{Myr}$) see the online extension of Figure 9 (top). While the majority of sink particle envelopes are reproduced well by Ulrich envelopes without singularities, this model under-predicts the density for some of the accreting objects. However, we cannot arbitrarily scale these particular profiles, because the mass M_* of the stellar object needs to drop for increasing ${\rho }_{u}^{\mathrm{env}}$ and constant infall rate ${\dot{M}}_{\mathrm{env}}$ (compare to Equation (23)).

Zoom In Zoom Out Reset image size

View extended figure (3 panels)

The power-law method e3, which fits the envelope profile, follows the slope of the radial profile of the outer envelope resolved by the SPH simulation, since it is an extrapolation of the gravitational collapse modeled in the D14 SPH simulations. We fit the slope of the profile γ and the density ${\rho }_{0}$ at the radius r₀, which we set as the distance of the closest SPH particle. To extrapolate (see Press et al. 1992) the profile of the density for different cells j within ${r}_{j}\lt {r}_{0}$, we are only using s3 circumstellar sites of the outer envelope (${r}_{j}\gt {r}_{0}$), which belong to the sink particles we want to fit the profile for. We only use sites from the outer envelope (see blue sample in Figure 9), since the inner density values (black) represent the background density of the outer envelope created by the SPH field close to the sink particle. In the bottom panel of Figure 9, we show the power-law fits of the accreting sink particles 11640, 110753, and 197247 in time-step 122 (6.109 $\mathrm{Myr}$), and we show the same for all accreting sink particles of this time-step in the online extension of Figure 9 (bottom). We limit the slope of the profile to the range $\gamma \in [0,-2]$, because for $\gamma =0$ the density follows the flat background density at the s3 circumstellar sites calculated in Section 3.3. The limit $\gamma =-2$ was chosen because it corresponds to the radial profile of the "singular isothermal sphere" in gravitational collapse (see Shu 1977). From Figure 9, we can see that the accretion stream (blue) is not spherically symmetric when the accreting sink particles sit in a close cluster of sink particles (11640, 197247) rather than being isolated (110753). We use a ${\chi }^{2}$ value (see Press et al. 1992), to quantify the quality of the envelope profile fit:

$$\begin{eqnarray}&&{\chi }^{2}=\underset{j=1}{\overset{{N}_{j}}{\displaystyle \sum }}\displaystyle \frac{{({\mathrm{log}}_{10}\rho ({r}_{j})-{\mathrm{log}}_{10}{\rho }_{\mathrm{fit}}({r}_{j}))}^{2}}{| \,{\mathrm{log}}_{10}\,{\rho }_{\mathrm{fit}}({r}_{j})| }.\end{eqnarray} \tag{ 29 }$$

We need N_j, which are all the sites that belong to the refined circumstellar sites of a sink particle with ${r}_{j}\gt {r}_{0}$, and ${\rho }_{\mathrm{fit}}({r}_{j})$, which is Equation (28) evaluated for the fitted parameters ${\gamma }_{\mathrm{fit}}$, and ${\rho }_{0,\mathrm{fit}}$ to calculate

$$\begin{eqnarray}&&\,{\mathrm{log}}_{10}\,{\rho }_{\mathrm{fit}}=\,{\mathrm{log}}_{10}\,{\rho }_{0,\mathrm{fit}}+{\gamma }_{\mathrm{fit}}\,{\mathrm{log}}_{10}({r}_{j}/{r}_{0}).\end{eqnarray} \tag{ 30 }$$

Spherically symmetric profiles produce a better fit and have lower ${\chi }^{2}$ values, as can be seen in the bottom panel of Figure 9. However, when comparing with method e2, which uses the rotational flattened envelope profile without singularity, we see that in the central region, the power-law density is several orders of magnitude higher.

While the Ulrich profile is consistent with the physical picture of a gravitational collapse with conservation of angular momentum, the gravitational collapse description in the D14 SPH simulations should reproduce power-law profiles on small scales. The power-law profiles in the D14 SPH simulations are limited by the isothermal accretion flow of $\rho \sim {r}^{-2}$ at larger radii ($r\gt 0.1\,\,\mathrm{pc}$). Additionally, within the accelerating collapsing shell ($r\lt 0.1\,\,\mathrm{pc}$), the velocity is not constant and the density profile is shallower, with $\rho \sim {r}^{-3/2}$ (Ulrich 1976; Shu 1977). This explains the offset of the radial density profile of the outer envelope from the simulations, compared to the rotationally flattened envelope profile and why the power-law profiles produce a better fit. Although the power-law profiles are consistent with the gravitational collapse description used in the D14 SPH simulations, the rotational flattened envelope profiles are consistent with the physical picture of small-scale gravitational collapse.

To explore the effects of the choice of refinement techniques on the synthetic observations, we will provide the results for both methods e2 and e3 in this paper.

We calculate the distance r_js, envelope density ${\rho }_{{js}}^{\mathrm{env}}$ of every s3 circumstellar site j and every accreting sink particle s for which we have already calculated the "background" density of method p3 ${\rho }_{{\rm{p}}3,j}$ in Section 3.3. We define that an envelope of a sink particle only contributes when ${r}_{{js}}\lt {r}_{\max }$, where r_max is the maximum distribution radius of s3 circumstellar sites that belong to a sink particle as described in Section 3.2, otherwise we set all elements of ${\rho }_{{js}}^{\mathrm{env}}=0$. We weigh (Press et al. 1992) the density ${\rho }_{{js}}^{\mathrm{env}}$ by

$$\begin{eqnarray}{w}_{j,s}=\left\{\begin{array}{ccc}1 & \forall & {r}_{{js}}\lt {r}_{0,s}\\ 1-\frac{{r}_{{js}}-{r}_{0,s}}{{r}_{\max }-{r}_{0,s}} & \forall & {r}_{0,s}\leqslant {r}_{{js}}\leqslant {r}_{\max },\\ 0 & \forall & {r}_{{js}}\gt {r}_{\max }\end{array}\right.\end{eqnarray} \tag{ 31 }$$

where ${r}_{0,s}$ is the radius of the closest sink particle, as described in Section 5.1. The "background" density ${\rho }_{{\rm{p}}3,j}$ is weighted with the counter weight

$$\begin{eqnarray}&&{\mathop{w}\limits^{\sim }}_{j,s}=1-{w}_{j,s}.\end{eqnarray} \tag{ 32 }$$

We set all weights ${w}_{j,s}$ and ${\mathop{w}\limits^{\sim }}_{j,s}$ to zero when the density ${\rho }_{{js}}^{\mathrm{env}}=0$.

For the Ulrich description of method e2, we sum the weighted density function to get the superposition envelope density ${\mathop{\hat{\rho }}}_{j,{\rm{p}}3}$ at every site j with

$$\begin{eqnarray}&&{\mathop{\hat{\rho }}}_{j,{\rm{p}}3}^{{\rm{e}}2}=\underset{s=1}{\overset{{N}_{s}}{\displaystyle \sum }}[{w}_{j,s}{\rho }_{{js}}^{\mathrm{env}}+{\rho }_{{\rm{p}}3,j}{\mathop{w}\limits^{\sim }}_{j,s}],\end{eqnarray} \tag{ 33 }$$

by using the normalization N_s, which is the number of accreting sink particles with envelopes. The summation is necessary, since the Ulrich envelope was constructed from the physical parameters, such as the infall rate ${\dot{M}}_{\mathrm{env}}$.

As explained above, for the power-law description e3, we take the maximum of the weighted density function when superimposing many envelopes

$$\begin{eqnarray}&&{\mathop{\hat{\rho }}}_{j,{\rm{p}}3}^{{\rm{e}}3}=\underset{s=1}{\overset{{N}_{s}}{\max }}[{w}_{j,s}{\rho }_{{js}}^{\mathrm{env}}+{\rho }_{{\rm{p}}3,j}{\mathop{w}\limits^{\sim }}_{j,s}].\end{eqnarray} \tag{ 34 }$$

When the superposition density is lower than the background density ${\rho }_{{\rm{p}}3,j}$, for example, in some cases the outer regions of the Ulrich envelope fits are too low, we set ${\mathop{\hat{\rho }}}_{j,{\rm{p}}3}={\rho }_{{\rm{p}}3,j}$. In Figure 10, we show the inward-extrapolated densities of the envelope profiles including superposition. For all objects in time-step 122 (6.109 $\mathrm{Myr}$), see the online extension of Figure 10. For the isolated object 110753, we only see the density extrapolation inwards, while for the objects 11640 and 197247, we see a superposition. For instance, in the upper panel in Figure 10, where we show the Ulrich description of method e2, we can see that for the non-isolated objects 11640 and 197247, the new profile is higher than the initial fit (green) with values approaching the initial fit at small radii.

Zoom In Zoom Out Reset image size

View extended figure (3 panels)

For the power-law envelope of method e3 in the lower panel of Figure 10, this effect was suppressed due to the fact that we take the maximum of the envelope densities rather than the sum, as shown in Equation (34). We can see that there is a smooth transition to other envelopes (black particles in objects 11640 and 197247) and also a seamless transition²² to the radial profile of other envelopes resolved by the SPH simulation (blue particles).

In Figure 11, we present a 2D slice through the density structure in order to show the improvement in resolution for the accreting sink particle 110753. The left panel shows the density evaluated on the grid where each cell corresponds to an SPH particle (s1), the central panel shows the improvement once the circumstellar cells are added (s3), and finally the panel on the right shows an inward extension of the envelope including the superposition of envelope densities ${\mathop{\hat{\rho }}}_{j,{\rm{p}}3}$ for a rotationally flattened envelope (e2).

Zoom In Zoom Out Reset image size

The simplest approach (hereafter, CM1), is not to add any circumstellar material on top of the density given by the SPH simulation, but still include circumstellar sites to ensure that the temperature profile around the sources is resolved. The radiative transfer calculation is then set up as listed in Table 2.

We extend the first approach by using the densities from the rotationally flattened envelopes (e2), we evaluated in Section 5.1, and also account for their superposition (see Section 5.2) and hereafter call this circumstellar material method, method CM2. We used the superimposed density ${\mathop{\hat{\rho }}}_{j,{\rm{p}}3}$ for sites of accreting sink particles; however, we set the superimposed envelope density ${\mathop{\hat{\rho }}}_{j,{\rm{p}}3}$ of sites j within the disk radius ${r}_{j}\lt {r}_{\max }^{d}=500\,\,\mathrm{au}$ (as defined in Equation (48)) to zero and rather include an Analytical YSO Model in the center. For illustration, see Figure 14 (middle).

Zoom In Zoom Out Reset image size

For the inner $500\,\,\mathrm{au}$, we set up the Analytical YSO Model using the radiative transfer code Hyperion, which is adding the spectral energy distribution (SED) of the YSO (star, envelope, cavity, and disk).

With the higher density material present, we need to adjust the mass of the source in the radiative transfer by the mass of the added circumstellar material (envelope, cavity, disk) between r₀ and ${r}_{\min }\approx 0$. Therefore, we calculate the envelope mass ${M}_{\mathrm{env}}^{{\rm{e}}2}$ by numerically integrating over the envelope density profile. We use the bipolar cavity setup described in Section 5.3 and set the outer radius of the envelope in the cavity description (Equation (36)) to ${R}_{\max }^{\mathrm{env}}=500\,\,\mathrm{au}$. The cavity removes about 2% from the mass of the envelope within ${R}_{\max }^{\mathrm{env}}$. We scale the disk density ${\rho }_{\mathrm{disk}}$ thus that the disk mass ${M}_{\mathrm{disk}}$ (from Equation (43)) is 1% of the stellar mass M_*. The other parameters of the Analytical YSO Models are set up with the parameters given in Section 5.4. The new mass of the stellar particle is then calculated as follows:

$$\begin{eqnarray}&&{M}_{* }={M}_{\mathrm{sink}}-({M}_{\mathrm{env}}^{{\rm{e}}2}-{M}_{\mathrm{cavity}}^{{\rm{e}}2})-{M}_{\mathrm{disk}}.\end{eqnarray} \tag{ 52 }$$

The method CM2 sets up with the parameters listed in Table 2. We use this setup to precompute the SED of the Analytical YSO Model. Since the Analytical YSO Model is not spherically symmetric, we compute the models for different viewing angles and determine an angle-averaged SED that we then use in the radiative transfer calculation. Therefore, we average the computed SED over 12 isotropically spaced viewing angles.

Method CM3 is similar to CM2 but with power-law envelopes instead of rotationally flattened Ulrich envelopes. The envelope mass ${M}_{\mathrm{env}}^{{\rm{e}}3}$ can be, therefore, calculated analytically as the spherical volume integral of the density function from Equation (27) between r₀ and ${r}_{\min }\approx 0$.

$$\begin{eqnarray}&&{M}_{\mathrm{env}}^{{\rm{e}}3}={\int }_{0}^{2\pi }d\phi {\int }_{-\pi }^{\pi }d\theta {\int }_{{r}_{\min }}^{{r}_{0}}{dr}\ {r}^{2}\sin \theta \rho (r)\end{eqnarray} \tag{ 53 }$$

$$\begin{eqnarray}&&=4\pi {\rho }_{0}{\int }_{{r}_{\min }\approx 0}^{{r}_{0}}{dr}\ {\left(\displaystyle \frac{r}{{r}_{0}}\right)}^{\gamma }{r}^{2}\end{eqnarray} \tag{ 54 }$$

$$\begin{eqnarray}&&\approx \displaystyle \frac{4\pi }{\gamma +3}\ {\rho }_{0}\ {r}_{0}^{3}\end{eqnarray} \tag{ 55 }$$

In Table 2, we summarize the setup parameters of method CM3. Note that the cavity removes about 4% from the mass of the envelope within ${R}_{\max }^{\mathrm{env}}$.

Above we presented the circumstellar material setup versions CM1 (star alone), CM2 (star and e2 Analytical YSO Models), and CM3 (star and e3 Analytical YSO Models).

Method CM1 is the simplest approach when setting up a radiative transfer model from simulations, because it does not require an inward envelope extrapolation. For instance, Offner et al. (2012) used this approach without circumstellar refinement because the circumstellar material was too optically thick to contribute to the flux that they were interested in.

However, in this work, we will present a parameter study. From now on, we will use methods CM1, CM2, and CM3 to set up the circumstellar material in the radiative transfer setup and therefore explore the effects of circumstellar material on the MIR flux.

In Figure 13, we show the properties of the Milky Way non-PAH dust from Weingartner & Draine (2001) with the re-normalization of abundances from Draine (2003) with ${R}_{V}=5.5$ and ${b}_{c}=3.0$, where R_V is the ratio of the visual extinction to reddening magnitude, and b_c is the concentration of carbon atoms in the medium. These dust properties reproduce the fluxes of the far-infrared well but lacks emission in the near-infrared (NIR) and MIR up to 20 μm.

We use an approximate method for the treatment of the emission from dust grains and PAH molecules, following Robitaille et al. (2012 and B. Draine 2016, private communication). We separate the dust size distribution into three grain species: 80.63% big grains (>200 Å); 13.51% of a smaller dust species within 200 Å and 20 Å, called very small grains (vsg); and 5.86% of the PAH molecules, called ultra-small grains (<20 Å, usg).

The very small and ultra-small grains are transiently heated as mentioned above, while the big grains are in local thermodynamical equilibrium (LTE), similarly to the non-PAH dust. When comparing the non-PAH dust (black) with the big grains (red) in Figure 13, we can see that the opacities are not that different, while the ultra-small grains (blue) and very small grains (green) have pronounced features in the MIR. The strong emissivity values of the ultra-small grains below 20 μm causes the PAH emission typically observed toward high-mass star-forming regions.

The method DT1 combines the radiative transfer temperature ${T}_{\mathrm{RT}}^{\mathrm{dust}}$, calculated from the non-PAH and PAH dust with the hydrodynamical temperature ${T}_{{\rm{p}}3}$.

$$\begin{eqnarray}&&{T}_{\mathrm{RT}}^{\mathrm{tot}}\,:{T}_{\mathrm{RT}}^{\mathrm{dust}}\mathrm{and}\ {T}_{{\rm{p}}3}.\end{eqnarray} \tag{ 56 }$$

Note that for the PAH dust only the big grains can heat thermally, and therefore, we couple only big grains with the hydrodynamical temperature. The temperatures are not directly summed, since this would not be physical. In Hyperion, the temperature is actually represented by the specific energy rate $\mathop{\dot{\epsilon }}$ (changes in absorption and emission process), which under the assumption of LTE can be transferred into a corresponding temperature. Therefore, the temperatures can be combined by converting them to $\mathop{\dot{\epsilon }}(T)$ and adding those energies before they are converted back.

We will further explore the effects when coupling ${T}_{\mathrm{RT}}^{\mathrm{dust}}$ with an isothermal temperature ${T}_{\mathrm{iso}}=18\,{\rm{K}}$ for the two different dust types in method DT2. The constant temperature of 18 K, on average, matches the temperature in relatively empty patches in the Galactic plane (see Paper II for more details), and is therefore a good choice for an ambient temperature

$$\begin{eqnarray}&&{T}_{\mathrm{RT}}^{\mathrm{tot}}\,:{T}_{\mathrm{RT}}^{\mathrm{dust}}\,\mathrm{and}\,{T}_{\mathrm{iso}},\end{eqnarray} \tag{ 57 }$$

Note that for the PAH dust, the temperature is only significant for the big grains, which satisfies LTE. Therefore, we couple only big grains with the isothermal temperature.

In the method DT3, as a sanity check, we only run the radiative transfer without combining with an extra temperature.

$$\begin{eqnarray}&&{T}_{\mathrm{RT}}^{\mathrm{tot}}\,:{T}_{\mathrm{RT}}^{\mathrm{dust}}\end{eqnarray} \tag{ 58 }$$

We produce ideal synthetic observations of the methods DT1, DT2, and DT3 using non-PAH and PAH dust at 20 μm and present them in Figure 14, where we also show a background subtracted image at WISE  22 μm of the Soul Nebula (Wright et al. 2010; Koenig et al. 2012), a real star-forming region, in the right panel.

We now focus on the non-PAH dust in the upper panels of Figure 14. We note that the non-PAH dust produces much too bright emission²⁵ in comparison to the real region, when coupling with the hydrodynamical temperature (see method DT1). The coupling of isothermal and radiative transfer temperature (DT2) produces drastically too low brightnesses (note the colorbar scaling factor 10⁻¹³) as does the radiative transfer temperature alone (DT3).

We found that the PAH dust in the lower panels of Figure 14 can reproduce the typical surface brightness of the Galactic star-forming region better. While the coupling with the hydrodynamical temperatures (see method DT1) causes even higher brightnesses than its non-PAH dust counterpart, neglecting the hydrodynamical temperatures produces better results: when combining the radiative transfer temperature with an ambient isothermal temperature (DT2), the brightness of the object is consistent with real observed regions. Also, method DT3 produces a comparable result to the real star-forming region but lacks emission in the outer regions for longer wavelengths.

From Figure 14, we can see that when combining the radiative transfer temperature with the hydrodynamical temperature (see method DT1), the brightness in the MIR is overestimated for non-PAH and PAH dust. The high brightness components are produced by the high hydrodynamical temperatures ${T}_{{\rm{p}}3}$ in the outer regions, which was coupled with the dust radiative transfer temperature in Equation (56). The gas temperatures are highest in the low-density regions, due to the Larson (2005) EOS (see Figure 2 and Section 2.2), but this is not observed for real star-forming regions. Although the temperatures of neutral SPH particles lie below the dust sublimation temperature of 1600 K (see Section 5.4), the brightness due to the hydrodynamical temperature is several orders too high, due to the erroneous assumption that the gas and dust are thermally coupled everywhere. In fact, the heating of the gas at low densities due to the Larson (2005) EOS used in the D14 SPH simulations, is physically due to the decoupling of the dust from the gas and therefore it does not make sense to use these high temperatures for the dust. Therefore, having explored the different ways of setting up the radiative transfer temperatures, we recommend using method DT2 by coupling the radiative transfer dust temperature ${T}_{\mathrm{RT}}^{\mathrm{dust}}$ with the ambient isothermal temperature ${T}_{\mathrm{iso}}$ for big grains of the PAH dust, which recovers the brightness features in the MIR much better.