A modular set of synthetic spectral energy distributions for young stellar objects

TL;DR

This paper introduces a new, flexible set of synthetic spectral energy distributions for young stellar objects, designed for initial modeling and broad parameter exploration across various evolutionary stages.

Contribution

The paper presents a new set of simplified, wide-parameter-space SED models for YSOs, improving coverage and suitability for far-infrared and sub-mm observations, with publicly available tools.

Findings

Models cover a wider range of YSO stages.

Models are less dependent on complex parameters.

Fitting examples demonstrate utility in real observations.

Abstract

In this paper, I present a new set of synthetic spectral energy distributions (SEDs) for young stellar objects (YSOs) spanning a wide range of evolutionary stages, from the youngest deeply embedded protostars to pre-main-sequence stars with few or no disks. These models include significant improvements on the previous generation of published models: in particular, the new models cover a much wider and more uniform region of parameter space, do not include highly model-dependent parameters, and include a number of improvements that make them more suited to modeling far-infrared and sub-mm observations of forming stars. Rather than all being part of a single monolithic set of models, the new models are split up into sets of varying complexity. The aim of the new set of models is not to provide the most physically realistic models for young stars, but rather to provide deliberately simplified models for initial modeling, which allows a wide range of parameter space to be explored. I present the design of the model set, and show examples of fitting these models to real observations to show how the new grid design can help us better understand what can be determined from limited unresolved observations. The models, as well as a Python-based fitting tool are publicly available to the community.

Tables

Table 1. Table 1: Parameter ranges sampled

Parameter	Symbol	Minimum	Maximum	Sampling
Stellar radius	$R_{\star }$	0.1 $R_{\odot }$	100 $R_{\odot }$	Log
Stellar temperature	$T_{\star }$	2000 K	30000 K	Log
Disk mass [dust]	$M_{\mathrm{disk}}$	${10}^{-8}$ $M_{\odot }$	$0.1$ $M_{\odot }$	Log
Disk inner radius	$R_{\min }^{\mathrm{disk}}$	$R_{\mathrm{sub}}$	1000 $R_{\mathrm{sub}}$	Log
Disk outer radius	$R_{\max }^{\mathrm{disk}}$	50 AU	5000 AU	Log
Disk flaring power	$\beta$	1	1.3	Linear
Disk surface density power	$p$	$-2$	$0$	Linear
Disk scaleheight	$h_{100AU}$	1 AU	20 AU	Log
Envelope density [dust]	${\rho }_0^{\mathrm{env}}$	${10}^{-24}$ g/cm3	${10}^{-16}$ g/cm3	Log
Envelope density power	$\gamma$	$-2$	$-1$	Linear
Envelope centrifugal radius	$R_c$	50 AU	5000 AU	Log
Cavity density [dust]	${\rho }_0^{\mathrm{cav}}$	${10}^{-23}$ g/cm3	${10}^{-20}$ g/cm3	Log
Cavity opening angle	${\theta }_0$	0∘	60∘	Linear
Cavity power	$c$	1	2	Linear

Table 2. Table 2: Model sets presented in this paper

Model set	Icon	Star	Disk	Envelope	Cavity	Ambient	Inner radius	Variables	Models
s---s-i		yes	…	…	…	…	…	2	10,000
sp--s-i		yes	passive	…	…	…	$R_{\mathrm{sub}}$	7	10,000
sp--h-i		yes	passive	…	…	…	variable	8	10,000
s---smi		yes	…	…	…	yes	$R_{\mathrm{sub}}$	2	10,000
sp--smi		yes	passive	…	…	yes	$R_{\mathrm{sub}}$	7	10,000
sp--hmi		yes	passive	…	…	yes	variable	8	10,000
s-p-smi		yes	…	power-law	…	yes	$R_{\mathrm{sub}}$	4	10,000
s-p-hmi		yes	…	power-law	…	yes	variable	5	10,000
s-pbsmi		yes	…	power-law	yes	yes	$R_{\mathrm{sub}}$	7	10,000
s-pbhmi		yes	…	power-law	yes	yes	variable	8	10,000
s-u-smi		yes	…	Ulrich	…	yes	$R_{\mathrm{sub}}$	4	10,000
s-u-hmi		yes	…	Ulrich	…	yes	variable	5	10,000
s-ubsmi		yes	…	Ulrich	yes	yes	$R_{\mathrm{sub}}$	7	10,000
s-ubhmi		yes	…	Ulrich	yes	yes	variable	8	10,000
spu-smi		yes	passive	Ulrich	…	yes	$R_{\mathrm{sub}}$	8	10,000
spu-hmi		yes	passive	Ulrich	…	yes	variable	9	10,000
spubsmi		yes	passive	Ulrich	yes	yes	$R_{\mathrm{sub}}$	11	40,000
spubhmi		yes	passive	Ulrich	yes	yes	variable	12	80,000

Equations

$$\rho (R,z,\varphi )={\rho }_0^{\mathrm{disk}}{\left(\frac{R_0}{R}\right)}^{\beta -p}\exp \left[-\frac{1}{2}{\left(\frac{z}{h(R)}\right)}^2\right],$$

$$h(R)=h_0{\left(\frac{R}{R_0}\right)}^{\beta }.$$

$$\rho (r)={\rho }_0^{\mathrm{env}}{\left(\frac{r}{r_0}\right)}^{\gamma },$$

$$\rho (r,\theta )={\rho }_0^{\mathrm{env}}{\left(\frac{r}{R_{\mathrm{c}}}\right)}^{-3/2}{\left(1+\frac{\mu }{{\mu }_0}\right)}^{-1/2}{\left(\frac{\mu }{{\mu }_0}+\frac{2{\mu }_0^2{R}_{\mathrm{c}}}{r}\right)}^{-1},$$

$${\rho }_0^{\mathrm{env}}=\frac{{\dot{M}}_{\mathrm{env}}}{4\pi {\left(G{M}_{\star }{R}_{\mathrm{c}}^3\right)}^{1/2}}.$$

$${\mu }_0^3+{\mu }_0\left(\frac{r}{R_{\mathrm{c}}}-1\right)-\mu \left(\frac{r}{R_{\mathrm{c}}}\right)=0.$$

$$z(R)=z_0{\left(\frac{R}{R_0}\right)}^c.$$

$$z_0$$

$$r_0$$

$$n_{\mathrm{models}}=10,000\times 2^{\max [n_{\mathrm{var}}-9,0]}.$$

$$p_{\mathrm{lin}}=\sqrt{\frac{Q^2+U^2}{I^2}},$$

$$p_{\mathrm{circ}}=\left|\frac{V}{I}\right|.$$

$$K\equiv \frac{P(D|M_1)}{P(D|M_2)}.$$

$$P(D|M)=\int p(\text{@vec}\theta |M)p(D|\text{@vec}\theta ,M)d\text{@vec}\theta .$$

$$P(D|M)=\sum _{j=1}^{N}P(\text{@vec}{\theta }_j|M)P(D|\text{@vec}{\theta }_j,M),$$

$$P(D|M)=\frac{1}{N}\sum _{j=1}^{N}P(D|\text{@vec}{\theta }_j,M).$$

$$P(D|\text{@vec}{\theta }_j,M)\propto \exp \left(-\frac{{\chi }^2}{2}\right).$$

$${\chi }^2-{\chi }_{\mathrm{best}}^2<3n_{\mathrm{data}},$$

$$P(D|\@vec{\theta_{j}},M)\propto\left\{\begin{array}[]{ccc}1&\rm{for}&\chi^{2}-\chi^{2}_{\rm best}\leq 3\,n_{\rm data}\\ 0&\rm{for}&\chi^{2}-\chi^{2}_{\rm best}>3\,n_{\rm data}\end{array}\right..$$

$${\chi }^2-{\chi }_{\mathrm{best}}^2\le 9n_{\mathrm{data}},$$

$$P(D|M)\propto \frac{N^{\mathrm{good}}}{N}.$$

Full text

Aims.

Methods.

*Results. Conclusions. *Astronomical databases: miscellaneous – Radiative transfer – Stars: formation – Stars: protostars

1 Introduction

Over the last two decades, increasingly sensitive and detailed surveys of our Galaxy at infrared wavelengths have resulted in a dramatic increase in the number of known young stellar objects (YSOs). The Spitzer/GLIMPSE survey at mid-infrared wavelengths contains tens of thousands of YSOs, with the sample of around 10,000 YSO candidates from robitaille:08:2413 representing only the brightest ones. Spitzer has also revealed thousands of YSOs in a number of individual regions, including over 8,000 YSO candidates in the survey of Cygnus X North (beerer:10:679), over 1,000 YSO candidates in the c2d survey of nearby star-forming regions (evans:09:321), and almost 3,500 YSO in the Orion A and B molecular clouds (megeath:12:192). The WISE survey (wright:10:1868), while less sensitive than Spitzer, covers the whole sky and has already unveiled thousands of YSOs (koenig:14:131). Herschel has both provided long-wavelength data for known YSO candidates, and has also revealed new extremely embedded protostars (e.g., stutz:13:36).

While the Hubble Space Telescope (HST) and the Atacama Large Millimeter Array (ALMA) allow the study of a smaller sample of YSOs with exquisite resolution and sensitivity, the vast majority of YSOs in the Galaxy will remain unresolved for a while yet. Nevertheless, with such a large sample of known YSOs and YSO candidates, we can make significant progress in quantifying the formation of stars across a wide range of environments. Extracting physical properties from limited – often unresolved – observations however requires radiative transfer modeling.

In robitaille:06:256, we presented a set of approximately 20,000 radiative transfer models, each of which had spectral energy distributions (SEDs) computed for ten viewing angles and 50 apertures. The models were set up by first sampling the central source mass between 0.1 and 50$M_{\odot}$, and the central source age between $10^{3}$ and $10^{7}$ yr (excluding post-main-sequence objects), and using evolutionary tracks to derive the source temperature and radius in each case. The parameters for the circumstellar environment, consisting of an accretion disk, infalling envelope, and bipolar cavities, were then sampled from ranges that were functions of the stellar mass and/or age. In total, 14 parameters were used to define the models.

In robitaille:07:328, we presented a tool to rapidly fit model SEDs from R06 to observations, correctly taking into account the effects of extinction and apertures. The aim of the tool was not to identify the best-fitting model, but rather the ensemble of models that could provide a good fit, thereby identifying likely ranges in parameters space.

The R06 models and the R07 fitting tool have been extensively used to model the SEDs of thousands of sources across many regions of star formation. These models have been invaluable in a number of studies to learn about the physical properties of young stars and regions. For example, forbrich:10:1453 carefully modeled the SEDs of a number of protostars in the IC348 and NGC2264 star-forming regions, demonstrating that it was crucial to include the intra-cluster extinction in order to reliably determine evolutionary stages of young sources in embedded clusters. Another example is the study by mottram:11:A149, which made use of the R06 models in order to determine bolometric fluxes and luminosities for over a thousand embedded massive stars in the Red MSX Source (RMS) survey.

In this paper, I present a new set of model SEDs for YSOs that significantly improves on the R06 models. Despite the success of the R06 models, a number of issues with those models remain, and these are described in Section 2. I then present an overview of the setup and methods for the new set of models in Sections 3 and 4 respectively, and present some initial results derived from the new models in Section 5. I discuss remaining caveats relating to the models in Section 6, and limitations inherent to SED modeling in Section 7.

2 Limitations of the previous models

In this section, I review the main limitations of the R06 models in order to justify the choices made for the new set of models. For a discussion of more general issues inherent in any SED modeling, the reader is referred to Section 7.

2.1 Model-dependent parameters

Of the 14 parameters defining the R06 models, some required specific assumptions to be made in order to be used in the analytical description of the model. For instance, the stellar mass and age had no direct impact on the SED, but instead, evolutionary tracks were needed to transform these into the stellar temperature and radius, which are the properties that actually have an effect on the SED. When fitting these models to observed SEDs, the parameters we can directly determine in a best-case scenario are the stellar temperature and luminosity (and therefore radius). The stellar mass and age are then only determined with the assumption of the evolutionary tracks. Similarly, the envelope infall rate was converted to an envelope density structure via the assumption of the Ulrich (1976) collapse model. Therefore, when one derives stellar masses, ages, and envelope infall rates from the R06 models, one is implicitly assuming a specific set of evolutionary tracks and a collapse model. However, users of the models may not necessarily realize the inherent assumptions, and most importantly cannot easily change these assumptions since they are built into the models.

2.2 Unevenly sampled parameter space

The R06 models are defined by sampling circumstellar environment parameters from ranges that are functions of the stellar mass and age. For example, the envelope infall rate is sampled from ranges that decrease over time, while the cavity opening angle is sampled from ranges that increase with time. The motivation for these was to restrict the parameter space coverage to regions that were thought to be realistic in order to minimize the computing time required, but the downside was that this led to correlations between parameters in the models before they were even applied to observations.

Examples include a correlation between the envelope infall rate $\dot{M}_{\rm env}$ and stellar mass $M_{\star}$ (since the envelope infall rate is sampled from a constant range in $\dot{M}_{\rm env}/M_{\star}$), between the envelope infall rate and the disk accretion rate (since both are sampled from ranges that decrease for larger ages), and an anti-correlation between infall rate and cavity opening angle (since, as mentioned above, the latter is sampled from a range that increases with time). These correlations can be seen in Figure 1 of robitaille:08:290. These built-in correlations mean that it is very difficult to look for such trends in samples of objects, since the built-in trends would mask any real trend.

2.3 Fixed model complexity

The R06 models, defined by 14 parameters, are often used without testing simpler models. For example, all of the models with envelopes in the R06 set also have disks, and furthermore the higher the envelope density, the higher the minimum disk mass found in the models. This means that there are, in fact, no deeply embedded protostellar models that do not have a disk. This in turn implies that it is not possible to use the models to try and find evidence for disks in embedded YSOs because there are no models to test the hypothesis that there is no disk.

2.4 Lack of cold dust at long wavelengths

The outer radii of the envelopes in the R06 models were defined as the radius where the temperature would drop to 30K if the envelope was optically thin to radiation. While for more embedded YSOs, the temperature may drop below 30K at the outer radius, those models only guarantee that they include all dust hotter than 30K. This in turn means that the models may not include enough cold material compared to observed sources, and may therefore show a deficit of long-wavelength emission.

The reason for this limitation was that the models were primarily designed with the Spitzer IRAC/MIPS and shorter wavelength data in mind, and at the time there were no surveys equivalent to the Herschel data, which now provides high-resolution long-wavelength fluxes for thousands of YSOs.

As a result, investigations using Herschel data have predictably found that the models are not always able to fit at the longest wavelengths. For example, Sewilo:10:L73 found that some of the Herschel sources in the Large Magellanic Cloud were not well fit by the R06 models beyond 100 $\mu$m. Because the R06 models do not include much dust below 30K, none of the model SEDs peak at wavelengths longer than 100 $\mu$m, but the coldest protostars observed with Herschel typically peak around 200 or 300 $\mu$m. While the R06 models were well-suited to Spitzer observations, new models are required now that Herschel observations of YSOs are common.

2.5 Poor signal-to-noise at long wavelengths

The R06 models were computed using the Monte-Carlo radiative transfer code developed by whitney:03:1079; whitney:03:1049; whitney:04:1177, which at the time computed the long-wavelength SED by sampling equal energy photon packets in the same way as for shorter wavelengths. However, traditional Monte-Carlo sampling done in this way results in the number of photons emitted being highest where the SED is brightest, and conversely lower where the SED is fainter, such as in the far-infrared and at mm wavelengths. As a result, most SEDs in the R06 set had poor signal-to-noise (S/N) at 1 mm and beyond, and some of the less embedded models even had poor S/N beyond 100 $\mu$m. While in many cases this has not been too problematic, in the sense that observations also have poor S/N where for the faintest sources, we ideally need models with high S/N at all wavelengths to make the most efficient use of available observations.

3 Model components and set-up

In this section, I present an overview of the new set of models, which addresses the issues described in Section 2, and includes a number of further improvements.

3.1 Design philosophy

The new models were not computed in a single monolithic set, but rather as several sets of models with increasing complexity. For instance, one of the sets consists of models with only a star with a surrounding disk, another set includes a disk and an envelope, but no bipolar cavities, and yet another set includes a disk, envelope, and bipolar cavities.

This modularity allows us to ask which model offers the best representation of the data, before even looking at the actual values of the parameters (see an example in Section 5.4). For example, a source might be fit by a complex model with a central source, a disk, and an envelope with bipolar cavities, but it might be equally well fit by a model with only a disk, indicating that there is no strong evidence one way or another for the presence of an envelope. It is important to assess not only the goodness of fit but also the simplicity or complexity of the model, since it is easy but not always meaningful to fit any set of data with an arbitrarily complex model. This design can also allow the available sets of models to grow over time.

The components used to generate the model sets are described in Section 3.2. For components containing free parameters, the free parameters were uniformly randomly sampled between a minimum and maximum value. The ranges of values used for each parameter are given in Table 1. By randomly sampling in uniform ranges, we can avoid correlations between parameters which were present in the R06 models. On the other hand, some of the combinations of parameters may be unphysical – since the definition of what might be considered physical will change over time (for example with stellar evolutionary tracks) the models presented here span a broad parameter space, and it is left to the user of the models to decide which models to ignore as being unphysical.

In general, models can be divided into two main categories. The first are models that aim to be as realistic as possible – for example in the case of radiative transfer, models with a realistic 3-d distribution of dust and with the best available dust model (with variable dust properties depending on location and temperature). The second category of models are models that are simpler but aim to provide insight into the effect of various physical processes, components, and so on. The collection of models presented in this paper is firmly in the second category: the aim is not to provide the most realistic model, but rather simple models that can be used to explore large regions of parameter space. This influences some of the decisions outlined in the following sections. I encourage users to use these models as a starting point for modeling observational data of young stars, but to follow this up with more detailed modeling if the observations cannot be reproduced with the simple models, or if spatially resolved data or spectra are available.

3.2 Components

As mentioned in §3.1 the new models consist of a combination of components which I describe in the following sections.

3.2.1 Central source

The central source, present in all models, was set to be a spherical source. Unlike the R06 models, the central source was not assigned a mass, and the stellar properties were not derived from evolutionary tracks. Instead, the central source was defined directly using a stellar radius $R_{\star}$ and effective temperature $T_{\star}$. While this does result in some of the models having unphysical combinations of $R_{\star}$ and $T_{\star}$, it allows the models to be independent of stellar evolutionary tracks. The idea is that users of the models can then – if needed – assume a specific set of tracks, decide which models are physical according to those tracks, and assign masses and ages to the stars.

The effective temperature of the source was used to select appropriate stellar photosphere models. For temperatures including and above 4,000 K, the photosphere models from castelli:04 were used, while for temperatures below 4,000 K, models computed with the PHOENIX code (brott:05) and intended for the GAIA mission111http://www.hs.uni-hamburg.de/EN/For/ThA/phoenix/gaia_info.html were used instead. Since the central source is not defined in terms of mass, we cannot calculate a surface gravity $\log{[g]}$. However, $\log{[g]}$ does not have a large impact on stellar photospheres in the range 3,000 K to 20,000 K, while above and below these temperatures, the largest difference depending on the choice of stellar photosphere model is generally of the order of 10% relative to the $\log{[g]}=4.0$ models. Therefore, the $\log{[g]}=4.0$ models were always used, with the caveat that these may be wrong by up to 10% for low and high temperatures, which affects mostly the near-infrared fluxes for non-embedded models.

3.2.2 Disk

In this set of models, only passive disks are included. For embedded YSOs, the effect of not explicitly including accretion is likely to be minimal because the heating from viscous dissipation in the disk is typically not very important, and the biggest effect is the increase in luminosity from the star. Since the star is embedded, the radiation from the star gets reprocessed and therefore the shape of the stellar spectrum is not important.

When modeling disks at near-infrared and longer wavelengths, the passive disks can still be used to model accreting disks for a similar reason: while we expect a little extra heating in the disk from the viscous dissipation, most of the accretion luminosity comes from the stellar surface, and again the increase in luminosity of the central source is the most important effect, so a model of a passive disk with a higher central luminosity star will likely be adequate in most cases.

Of course, the lack of accretion does mean that UV and optical fluxes for non-embedded sources with strong accretion cannot reliably be modeled, since none of the passive disk models will produce the adequate excess UV and optical emission typically observed toward accreting sources.

Disks in hydrostatic equilibrium are expected to be flared (shakura:73:337), and at earlier times, when the dust is well coupled to the gas, the dust in the disk follows the same structure. Over time, as the dust settles to the mid-plane, the effective flaring for the disk may change. The flaring and the scaleheight of the disk is therefore parametrized such that it covers the range of flaring from hydrostatic disks to flat disks.

The density distribution of the passive flared disks is given in cylindrical polar coordinates $(R,z,\phi)$ by

[TABLE]

where $\rho_{0}^{\rm disk}$ is defined by the disk dust mass $M_{\rm disk}$, $p$ is the surface density radial exponent, $\beta$ is the disk flaring power, and the disk scaleheight $h(R)$ is given by

[TABLE]

The disk is truncated at the inner radius $R^{\rm disk}_{\rm min}$ and the outer radius $R^{\rm disk}_{\rm max}$.

The free parameters varied for all models with a disk are the disk dust mass $M_{\rm disk}$, outer radius $R^{\rm disk}_{\rm max}$, flaring exponent $\beta$, surface density exponent $p$, and scaleheight $h_{0}$. For some of the models, the disk inner radius $R^{\rm disk}_{\rm min}$ was also varied (as described in §3.3).

3.2.3 Envelope

Two types of envelopes were included in the models – the first are spherically-symmetric power-law envelopes, and the second are rotationally flattened envelopes. The reason for including both types of envelopes is because this will allow users to investigate the constraints on the envelope structure from observations: while a model with a more complex rotationally flattened envelope may fit well, it is important to also check whether a model with a simpler spherically symmetric envelope can also fit. The spherically symmetric models also have more flexibility as to what power-law to use for the radial dependence of the density.

The power-law envelope density structure is given by

[TABLE]

where $\rho_{0}^{\rm env}$ is the density of the envelope at the radius $r_{0}$, and serves as the scaling for the envelope density, and $\gamma$ is the radial power of the density.