<figure><img src="image_1.png" /></figure>

# The Nature of Class I Sources: Periodic Variables in Orion

Elizabeth A. Ellithorpe$^{1,2}$, Gaspard Duchene$^{2,3}$✉, and Steven W. Stahler$^2$

Received 2019 July 8; revised 2019 September 17; accepted 2019 September 18; published 2019 October 31

## Abstract

We present a quantitative, empirically based argument that at least some Class I sources are low-mass, pre-main-sequence stars surrounded by spatially extended envelopes of dusty gas. The source luminosity arises principally from stellar gravitational contraction, as in optically visible pre-main-sequence stars that lack such envelopes. We base our argument on the fact that some Class I sources in Orion and other star-forming regions have been observed by *Spitzer* to be periodic variables in the mid-infrared, and with periods consistent with T Tauri rotation rates. Using a radiative transfer code, we construct a variety of dust envelopes surrounding rotating, spotted stars, to see whether an envelope that produces a Class I spectral energy distribution at least broadly matches the observed modulations in luminosity. Acceptable envelopes can be either spherical or flattened and may or may not have polar cavities. The key requirement is that they have a modest equatorial optical depth at the *Spitzer* waveband of 3.6 µm, typically $\tau_{3.6} \approx 0.6$. The total envelope mass, based on this limited study, is at most about 0.1 $M_\odot$, less than a typical stellar mass. Future studies should focus on the dynamics of the envelope, to determine whether material is actually falling onto the circumstellar disk.

*Unified Astronomy Thesaurus concepts: Circumstellar matter (241); Pre-main sequence stars (1290); Protostars (1302); Starspots (1572); T Tauri stars (1681); Herbig Ae/Be stars (723)*

*Supporting material: figure set, machine-readable table*

## 1. Historical Introduction

The evolutionary status of young stars embedded in dusty gas remains problematic. Traditionally, researchers have utilized the observed spectral energy distribution (SED) of a source to infer the age of its central star (e.g., Evans et al. 2009). The youngest sources, dubbed Class 0, lack both optical and near-infrared flux, and their SEDs continue to rise beyond 100 µm (André et al. 1993). More quantitatively, Class 0 sources have the lowest bolometric temperature, $T_{bol}$, where $T_{bol}$ is a global measure of the envelope's mean dust temperature (Myers & Ladd 1993). They also have submillimeter luminosities ($\lambda > 350$ µm) greater than about 0.5% of the bolometric value, indicating that the mass of the star's envelope exceeds that of the star itself (André et al. 1993).

Within this same empirical framework, the next-youngest category is Class I. Here, there is detectable near-infrared flux, and the SED rises from 2.2 to 20 µm. The envelope appears to contain less mass than the star, since the millimeter luminosity is a smaller fraction of $L_{bol}$ than in Class 0 sources. By now, it is common practice for observers to refer to both Class 0 and Class I objects as “protostars” (e.g., Furlan et al. 2016) and to assume that the latter are, in some sense, more evolved than the former.

Exactly what this sense is remains unclear. Over the past 15 yr, many observers have probed various features of Class I sources, in part to address this evolutionary issue. White & Hillenbrad (2004) observed 15 Class I stars in Taurus optically, via scattered light. Spectroscopy indicated that both the stars' effective temperature and rotation speeds were indistinguishable from the optically revealed Class II sources, long identified as more evolved, pre-main-sequence stars. However, the target sources could have been a relatively old subset of Class I sources, since they had envelopes that were porous

enough to allow the escape of optical photons. Using a broader sample, and focusing on the near-infrared spectrum, Covey et al. (2005) found a larger mean rotation speed than in Class II objects. On the other hand, Doppmann et al. (2005) found the stellar luminosity of these objects to be very similar to classical T Tauri stars. It should be noted that obtaining $L_*$ in any source pivots on a derived value for the extinction. Unfortunately, obtaining this extinction is highly problematic in Class I, a point emphasized by Beck (2007).

Prato et al. (2009) and Connelley & Greene (2010) estimated Class I mass accretion rates, using the observed flux in the near-infrared $Br\gamma$ line and an empirical relationship initially established for T Tauri stars (Muzerolle et al. 1998). Most of the resulting values for the accretion luminosity $L_{acc}$ were far less than the total luminosity of the source. But even this point remains controversial. To cite one specific example, Beck (2007) used $Br\gamma$ emission to estimate $L_{acc}/L_{bol} = 0.008$ for the source IRAS 04016+2610 in Taurus. Furlan et al. (2008), who fit a multi-parameter model to the star's SED, derived $L_{acc}/L_{bol} = 0.90$.

Modeling the SED yields, at least potentially, both the mass of the envelope and the rate at which gas is falling onto the star. Matching the observed SEDs with those from models of collapsing, rotating clouds suggests Class I infall rates of order $10^{-6} M_\odot \text{ yr}^{-1}$ (Whitney et al. 2003), two orders of magnitude higher than those typically obtained for Class II objects (Hartmann et al. 2016). On the other hand, the infall rates derived via infrared spectroscopy are lower by several orders of magnitude (e.g., Prato et al. 2009). There is still no generally accepted resolution to this discrepancy.

Fortunately, the temporal behavior of Class I sources provides an important clue to the nature of the central object and therefore to the evolutionary state. Low-mass, pre-main-sequence stars that do not reside in dusty envelopes are well

<footnote>$^1$Physics and Astronomy Department, University of Oklahoma, Norman, OK 73019, USA; ellithorpe@ou.edu</footnote>

<footnote>$^2$Astronomy Department, University of California, Berkeley, CA 94720, USA</footnote>

<footnote>$^3$Université Grenoble-Alpes, CNRS Institut de Planétologie et d'Astrophysique (IPAG), F-38000 Grenoble, France</footnote>
known for their variability, which can be either erratic or quasi-
periodic (e.g., Cody et al. 2014). In the latter case, hot or cool
spots on the stellar surface, rotating in and out of view,
plausibly explain the luminosity fluctuations, although asym-
metric disk emission can also contribute at infrared
wavelengths.

In recent years, *Spitzer* monitoring campaigns of star-
forming regions have revealed that some Class I sources also
display quasi-periodic photometric variations (e.g., Morales-
Calderón et al. 2011; Wolk et al. 2015), further suggesting that
these objects could be spotted pre-main sequence stars.
However, the interpretation of this finding is not straightfor-
ward, as the infrared light is reprocessed by the circumstellar
envelope. Thus, one might expect that the spot emission would
be redistributed in a more isotropic fashion and toward longer
wavelengths, suppressing the characteristic near-infrared
photometric signal. This effect strongly depends on the
envelope structure and total optical depth. Therefore, to
properly interpret the implications of the periodic variability
of Class I sources, it is necessary to assess how deeply
embedded they truly are.

The present paper develops this line of argument quantita-
tively. We employ a radiative transfer code to explore
numerical models of dust envelopes surrounding rotating,
spotted stars. We test these models against observational data,
consisting of a set of periodically varying Class I sources in
Orion. We judge models to be acceptable only if they yield
both Class I SEDs and infrared light curves that resemble those
of our observed sample. This double constraint allows us to
characterize important features of the envelopes' density
distribution in more detail than previous studies.

Our main result is that Class I envelopes, while sufficiently
opaque to block starlight from reaching the observer directly,
must have a modest mid-infrared optical thickness ($\tau_{3.6} \approx 0.6$)
to preserve an intrinsic periodic signal.

Flattened envelopes similar to those in rotating, infalling
clouds are acceptable, if the envelope mass is sufficiently low.
So are purely spherical clouds with power-law density profiles,
as empirically found in Class 0 envelopes (e.g., Tobin et al.
2015), again provided that the mid-infrared optical thickness is
modest. Neither SED analysis nor photometric light curves
allow us to constrain the detailed structure of the envelope,
which can only be determined via future high-resolution
imaging studies.

In any event, since the observed periodicity of the Orion
Class I sources is compatible with stellar rotation rates, our
findings add weight to the possibility that Class I sources
generally are embedded pre-main-sequence stars.

We begin, in Section 2 below, by exploring which dust
envelopes yield appropriate Class I SEDs. From a purely
practical standpoint, we require such an envelope for our
analysis of the periodic behavior of an embedded, spotted
Class I object. To elucidate the best-fitting model envelope, we
use a radiative transfer code to produce synthetic SEDs, which
we compare to an empirical, Class I template. In Section 3, we
review the *Spitzer* observations of periodically varying Class I
sources in Orion. We display the SEDs and construct light
curves for all stars. We account for these light curves by
assuming that there are single spots located on the rotating
stellar surface. We then utilize our previously constructed set of
envelopes yielding proper SEDs to see which models also
produce satisfactory light curves. We also review, for

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Luminosity-normalized SEDs generated by three TSC-type envelopes, along with the empirical Class I template SED. Shown are SEDs in the best-fitting case, along with those from the lowest and highest acceptable envelope masses. Parameters for all three models are listed in Table 1.</figcaption></figure>

comparison, the observations of periodic Class II and III
sources; our Orion Class I sample is more akin to the former
than to the latter. Section 4 summarizes our findings and
discusses the broader astrophysical implications of this work.

**2. Class I Envelopes**

2.1. *Modeling Strategy*

We first briefly examine which dusty envelopes surrounding
a young star produce an SED typical of Class I sources. To this
end, we compare synthetic SEDs to observed ones. Empiri-
cally, there is a rather broad range of SED morphologies that all
share the defining Class I property of a rising flux between 2
and 20 µm. We have therefore constructed a representative, or
template, SED by suitably averaging photometric data from 21
sources in Taurus, the star-forming region with the best-studied
population of Class I sources. To combine the SEDs, we first
normalized each to the bolometric luminosity of the system, so
as to focus on the intrinsic shape of the flux distribution. We
refer the reader to Appendix A for the remaining details of our
procedure for constructing the template. We also built a
template based on Orion Class I sources, which is similar to the
Taurus one, albeit with some biases introduced by crowding
and confusion. We therefore adopt the Taurus template,
illustrated in Figures 1 and 2, for our analysis.

For constructing models of Class I sources, we assume three
components: a pre-main sequence star at the center, a
circumstellar disk, and an extended envelope. The parameters
describing each component are presented in the remainder of
this section. Each synthetic SED was obtained using the 2D
version of the radiative transfer code MCFOST (Pinte et al.
2006, 2009). To evaluate whether a given model adequately
matches the template Class I SED, we used a standard χ²-test.
We account for the asymmetric error bars associated with the
template by adopting the upper or lower uncertainty depending
on whether the model over- or underpredicts the flux of the
template at a given wavelength. We consider models yielding
reduced χ²<sub>red</sub> ≤ 1 to be satisfactory matches to the template
SED. (Of the objects used in defining the template, two-thirds
pass this criterion when compared to the template itself.)

The parameter space describing a Class I system is
impractically large to perform a thorough exploration. Besides
<figure><img src="image_3.png" /><figcaption>**Figure 2.** SEDs generated by three power-law envelopes ($\gamma = 1.0$), along with the template Class I SED. Shown are SEDs from the optimal, lowest-mass, and highest-mass model envelopes.</figcaption></figure>

the sheer number of parameters being a limiting factor in how detailed our analysis can be, SEDs are notoriously degenerate, making convergence on a single model difficult. We therefore decided to fix as many parameters as possible and only allow the total disk and envelope masses to vary in our investigation of acceptable models of Class I systems. In all models, we assume that dust composes 1% of the gas mass. To evaluate the SED of the system, we also need to select a viewing angle. We choose a polar angle of $\theta = 57^{\circ}$, which is the mean value for an envelope with random orientation.

For the central star, we adopt a total luminosity of $0.9 L_{\odot}$. We further assume that the spectrum of the object is well represented by a 4000 K NextGen stellar atmospheric model (Allard et al. 2012), typical for Class II and III sources in star-forming regions. Thus, the stellar radius is $2 R_{\odot}$. For our radiative transfer calculation, we need not specify the stellar mass. However, we do require this quantity in one of our envelope models (see below), and for that purpose we choose the representative value $M_{*} = 0.5 M_{\odot}$.

It has long been known that, in regions as young as Orion, the vast majority of optically visible pre-main-sequence stars are surrounded by circumstellar disks (Haisch et al. 2001). We therefore assume that our more embedded stars also have such disks. Our model disks have a surface density that varies as $r^{-3/2}$, a 100:1 gas-to-dust mass ratio, and a standard flared structure. Thus, the scale height $H$ is given by $H_0 (r/r_0)^{6/5}$, where $H_0 = 10$ au and $r_0$ is the outer radius of 100 au. The inner radius is 0.073 au, which is the representative value adopted by Whitney et al. (2003) for the location where dust is thermally sublimated. The dust itself is characterized by a 3.5:1 silicate/graphite mixture and a power-law size distribution extending from 0.03 to 1 mm, with an index of $-3.5$ (Mathis et al. 1977). Finally, we assume that the disk itself is passively irradiated by the star, i.e., we do not include any additional luminosity either from infall onto the disk surfaces or from accretion within that structure.

In exploring the envelope structure, we consider two model families (see Table 1) as in some previous studies of Class I sources (e.g., Furlan et al. 2008). Both families have the same inner and outer radii, 0.073 and 5000 au, respectively. Here, we are guided by the SED study of Whitney et al. (2003). We further note that the outer radius has no significant influence on the resulting SED, whereas the inner radius must not exceed

$\sim 1$ au in order to match the Class I template SED in the 1–5 $\mu$m range. The dust in our model envelopes has the same compositional mixture as for the disk, but a size distribution extending from 0.03 to 1 $\mu$m with a $-4.5$ power law. We again assume a 100:1 gas-to-dust ratio. These parameters were chosen so as to produce an extinction law that matches observations in the interstellar medium and that was used in Class I models by Whitney et al. (2003).<fnref n="4" />

## 2.2. TSC Envelopes

For over 30 yr, many researchers seeking to model nonspherical envelopes around young stars have utilized the classic study of Cassen & Moosman (1981). These authors calculated the density distribution in the deep interior of a collapsing, rotating cloud. The cloud itself is the highly idealized singular, isothermal sphere, whose self-similar collapse was previously studied by Shu (1977). Cassen & Moosman (1981) found an analytic expression for the density as a function of radius $r$ and polar angle $\theta$:

$$ \rho(r, \theta) = \frac{\dot{M}_{\text{env}}}{4\pi(GM_{*}R_{c}^{3})^{1/2}} \left(\frac{r}{R_{c}}\right)^{-3/2} \times \left(1 + \frac{\mu}{\mu_{0}}\right)^{-1/2} \left(\frac{\mu}{\mu_{0}} + \frac{2\mu R_{c}}{r}\right)^{-1} \quad (1) $$

In this equation, the quantity $R_c$ represents the centrifugal radius, i.e., the outermost point where fluid trajectories land in the equatorial plane. The density is peaked toward the central star, varying with radius as $r^{-1/2}$ for $r \ll R_c$ and as $r^{-3/2}$ for $r \gg R_c$. A second quantity appearing in the equation is $\dot{M}_{\text{env}}$, the total mass accretion rate onto the disk and star, whose mass is $M_{*}$. The expression also contains the variables $\mu \equiv \cos \theta$ and $\mu_0 \equiv \cos \theta_0$, where $\theta_0$ is the polar angle at which each fluid element leaves the rarefaction wave. Since the element subsequently traces a parabola, $\mu$ is a function of both $r$ and the parameters $\mu_0$ and $R_c$. This functional dependence is given implicitly by

$$ \mu_0^3 + \mu_0 \left( \frac{r}{R_c} - 1 \right) - \mu \left( \frac{r}{R_c} \right) = 0. \quad (2) $$

The studies by Cassen & Moosman (1981) and Shu (1977) provide expressions for both $\dot{M}_{\text{env}}$ and $R_c$. The first is a function only of the parent cloud temperature (assumed to be uniform), while the second also depends on the cloud's rotation rate and the time since the start of collapse. In practice, researchers have treated $\dot{M}_{\text{env}}$ and $R_c$ as free parameters, thereby generating families of envelope solutions. These generalized models, which we designate as TSC, are not physical solutions to the problem of a rotating, collapsing cloud. They are simply a set of nonspherical envelopes useful for characterizing observed sources.

We follow this practice here, while noting, for example, that $\dot{M}_{\text{env}}$ cannot be interpreted literally as the infall rate in the envelope; it is simply a scale factor in the expression for envelope density. In any event, we implemented Equations (1) and (2) to construct a dust envelope model for use in MCFOST. Following Whitney et al. (2003), we set $R_c = 100$ au. With these values for $R_c$ and the inner and outer radii of the

<footnote n="4">⁴ More recent observations indicate that the extinction curve in molecular clouds is flatter than that in the interstellar medium for wavelengths $\lambda \gtrsim 3 \mu$m (Flaherty et al. 2007).</footnote>
envelope, the total mass of the envelope is proportional to the mass accretion rate in the envelope, with

$$M_{\text{env}} = 0.19 \left( \frac{\dot{M}_{\text{env}}}{5 \times 10^{-6} M_{\odot} \text{yr}^{-1}} \right) M_{\odot}. \quad (3)$$

The solid curve in Figure 1 displays the best-fit SED, along with the discrete points and error bars corresponding to our Class I template. For this model, the envelope mass was $M_{\text{env}} = 0.050 M_{\odot}$, and the disk mass was $M_{\text{d}} = 2 \times 10^{-3} M_{\odot}$. Note that $M_{\text{env}} < M_{\text{*}}$, a property that we found to hold generally, in line with previous studies of the SED of Class I sources (e.g., Whitney et al. 2003; Furlan et al. 2016). The optimal TSC model has $\chi^2_{\text{red}} = 0.23$.

A key property of the optimal model is its modest optical depth at mid-infrared wavelengths. Referring to the Spitzer waveband of 3.6 µm, we find $\tau_{3.6} = 0.8$ at our viewing angle of 57°. The corresponding figure is 0.6 (1.9) along the pole (equator). As we will see, limited optical depth turns out to be the essential property for reproducing a Class I-type SED.

To illustrate this point further, we also display in Figure 1 models that have the lowest and highest acceptable envelope masses; these both have $\chi^2_{\text{red}}$-values close to unity. The dashed SED ($\chi^2_{\text{red}} = 0.95$) was obtained from a model with $M_{\text{env}} = 0.027 M_{\odot}$ and $M_{\text{d}} = 3 \times 10^{-3} M_{\odot}$. Here, $\tau_{3.6} = 0.4$ at $\theta = 57°$. With such a relatively low optical depth, the short-wavelength radiation from the star and inner disk partially penetrates the envelope, lifting the SED above the template in this regime. If we were to include an outflow cavity, as was done by Whitney et al. (2003) in several models, then the short-wavelength SED would be even more elevated.

Finally, the dotted curve in Figure 1 corresponds to our highest acceptable envelope mass, $M_{\text{env}} = 0.1 M_{\odot}$. In this case, the SED has $\chi^2_{\text{red}} = 1.04$. The disk mass is actually lower than for the optimal model, $M_{\text{d}} = 5 \times 10^{-4} M_{\odot}$. The SED mimics the template quite well at short wavelengths but rises above it at longer wavelengths, representing absorption and reemission of radiation by the relatively optically thick envelope.

We see that dust envelopes constructed according to the TSC model yield SEDs whose shapes can match typical Class I sources. However, the luminosities predicted by these same models are too large, by an order of magnitude. Thus, the accretion luminosity corresponding to our optimal TSC model, $L_{\text{acc}} = GM_{*}\dot{M}_{\text{env}}/R_{*}$, is $10 L_{\odot}$, provided that we interpret the parameter $\dot{M}_{\text{env}}$ as an infall rate. In contrast, the median bolometric luminosity of the Taurus Class I sources used in generating the template is only $0.8 L_{\odot}$. This discrepancy is a manifestation of the venerable “luminosity problem” (Kenyon & Hartmann 1995). However, the paradox disappears if we recognize, as already stressed, that $\dot{M}_{\text{env}}$ in TSC models is simply a parameter governing the envelope mass and does not reflect the true infall rate, if any.

### 2.3. Power-law Envelopes

TSC models resemble, at least qualitatively, rotating clouds undergoing inside-out collapse. However, the inferred envelope masses that yield Class I SEDs are much smaller than observed values for dense cores. We noted that a TSC envelope yields an acceptable SED only when the mid-infrared optical thickness is modest. It is this feature, rather than the dynamics of the original TSC theory, that is essential.

We bolster this statement by considering envelopes that, unlike TSC models, have no polar variation in density. That is, we assume the density to be spherically symmetric and to vary as a power law in radius: $\rho \propto r^{-\gamma}$. In addition to their simplicity, such models have some empirical motivation. Resolved millimeter maps of Class 0 sources find density variations that are at least consistent with power laws. Furthermore, there is no indication of a change in the index at an $R_c$ value of order 100 au (e.g., Chiang et al. 2012; Tobin et al. 2015; Persson et al. 2016). Unfortunately, such finer mapping is not yet feasible for the lower-mass Class I envelopes, so the finding for less evolved objects is only suggestive.

Keeping the previous star and disk models, we now construct power-law envelopes that still have inner and outer radii of 0.07 and 5000 au, respectively. For the exponent $\gamma$, we choose 1. This value is broadly consistent with the millimeter observations (e.g., Tobin et al. 2015; Persson et al. 2016) and differs from the $\gamma$ of 1/2 that characterizes the deep interior of TSC models.

We display the SED generated by the optimal power-law envelope as the solid curve in Figure 2. It is evident visually that the match to the template SED is quite close; the associated $\chi^2_{\text{red}}$ value is 0.14, even smaller than for the optimal TSC model. In this case, $\tau_{3.6} = 1.5$ at all polar angles. The mass of the envelope is now $0.13 M_{\odot}$, while that of the disk is $1 \times 10^{-3} M_{\odot}$.

The dashed curve in Figure 2 shows the SED for the lowest-mass envelope that still yields an acceptable match to the Class I template, i.e., with a $\chi^2_{\text{red}}$ near unity. Here, $M_{\text{env}} = 0.05 M_{\odot}$, $M_{\text{d}} = 1 \times 10^{-3} M_{\odot}$, and $\chi^2_{\text{red}} = 0.94$. The optical thickness $\tau_{3.6}$ is 0.6. As in the analogous TSC model, the short-wavelength portion of the synthetic SED rises markedly above the template. In the maximum-mass envelope, whose SED is shown as the dotted curve, the short-wavelength portion is suppressed, while the SED rises in the mid-infrared region with a marked silicate absorption feature that is present in only some Class I sources (Furlan et al. 2008). Here, $M_{\text{env}} = 0.28 M_{\odot}$, $M_{\text{d}} = 3 \times 10^{-3} M_{\odot}$, and $\chi^2_{\text{red}} = 0.98$. The relevant optical depth is now $\tau_{3.6} = 3.2$.

To compare our model envelope masses, both TSC and power law, with observed ones, we note that the K-band extinction is a convenient proxy for the latter. The dust extinction law adopted for the envelope is such that $A_K/\tau_{3.6} \approx 3.3$. Thus, the typical line-of-sight extinction associated with Class I envelopes is a few magnitudes at K. This is fully consistent with the empirical estimates based on spectral and photometric analysis of known Class I systems (Beck 2007; Prato et al. 2009; Connelley & Greene 2010; Furlan et al. 2016).

## 3. Periodic Class I Sources in Orion

### 3.1. Source List

This study was initially motivated by the identification of several periodically variable Class I sources in Orion, based on a 40-day 3.6 and 4.5 µm photometric monitoring campaign with Spitzer (Morales-Calderón et al. 2011). To confirm the objects’ class and status as periodic variables, we first cross-referenced the sample of periodic sources in Morales-Calderón et al. (2011) with the SED classification database subsequently presented by Megeath et al. (2012). We thus identified nine
<caption>Table 1</caption>

Parameters of Synthetic Models That Assume a TSC Envelope (Top Part of the Table) and a Power-law Envelope (PL, with $\gamma = 1$; Bottom Part of the Table)

<table><thead><tr><th>Model</th><th><i>M</i><sub>env</sub> (<i>M</i><sub>☉</sub>)</th><th><i>M</i><sub>disk</sub> (<i>M</i><sub>☉</sub>)</th><th>τ<sup>env</sup><sub>3.6</sub></th><th>χ<sup>2</sup><sub>red</sub></th></tr></thead><tbody><tr><td>TSC low-mass</td><td>2.7 10<sup>-2</sup></td><td>3 10<sup>-3</sup></td><td>0.3–1.0</td><td>0.95</td></tr><tr><td>TSC optimal</td><td>5 10<sup>-2</sup></td><td>2 10<sup>-3</sup></td><td>0.6–1.9</td><td>0.23</td></tr><tr><td>TSC high-mass</td><td>1 10<sup>-1</sup></td><td>5 10<sup>-4</sup></td><td>1.1–3.8</td><td>1.04</td></tr><tr><td>PL low-mass</td><td>5 10<sup>-2</sup></td><td>3 10<sup>-3</sup></td><td>0.6</td><td>0.94</td></tr><tr><td>PL optimal</td><td>1.25 10<sup>-1</sup></td><td>1 10<sup>-3</sup></td><td>1.5</td><td>0.14</td></tr><tr><td>PL high-mass</td><td>2.75 10<sup>-1</sup></td><td>3 10<sup>-5</sup></td><td>3.2</td><td>0.98</td></tr></tbody></table>

**Note.** The total gas masses in the envelope and disk are listed in Columns (2) and (3), respectively. Column (4) indicates the 3.6 µm optical depth through the envelope; for the TSC envelope, the lowest (highest) optical depth is found at the pole (equator). Column (5) presents the $\chi^2_{red}$ when compared to the Class I template.

possible sources as potential Class I periodic variables. We subsequently discarded two of these sources, as their SEDs are clearly inconsistent with Class I sources. Instead, they are Class II or III sources that were misclassified in Megeath et al. (2012). The remaining seven objects, on which we henceforth focus, are listed in Table 2. To ensure that these seven objects are truly representative of Class I sources, we consulted the color–color diagrams of Megeath et al. (2012) and placed our objects in them using photometry data. The one diagram in which we can display all our objects is shown as Figure 3. The near- and mid-infrared colors of our target stars span the range of Class I sources as identified by Megeath et al. (2012). In the absence of systematic differences in color, we conclude that these sources are indeed representative of the larger population of Class I sources in Orion.

Our next step was to construct more complete SEDs, using both Two Micron All Sky Survey (2MASS) JHK photometry (Skrutskie et al. 2006) and the Spitzer photometry published in Megeath et al. (2012). Additional photometry data for sources 5 and 7 were obtained from the HOPS survey in Furlan et al. (2016). The SEDs of all seven sources are shown in Figure 4. In all cases, the SEDs are appropriately identified as Class I, since the flux density either steadily increases toward longer wavelengths or is flat across the near- to mid-infrared regime; this behavior is also consistent with the Class I template we have assembled (recall Figures 1 and 2).

We also inspected the 2MASS, Spitzer, and WISE (Wright et al. 2010) images of each target to search for potential source confusion. Specifically, we wanted to test the possibility that some of the Class I periodic sources are actually composed of two unrelated objects that lie in close proximity in the sky, a periodically variable star (possibly a Class II source in Orion) and an infrared-bright source. Given its relatively low angular resolution (several arcseconds, depending on wavelength), Spitzer would not be able to resolve the two objects and would thus erroneously assign the periodicity and infrared excess to the same object. All targets were detected in the 2MASS images, and none show any sign of confusion. It is still possible that some binary systems, such as a Class II source and a so-called “infrared companion” at close separation (Koresko et al. 1997), could still escape inspection in this manner. However, the low occurrence rate of such systems indicates that they are unlikely to be present in all cases. In addition, Sources 2, 6, and 7 have been assessed to be single in near-infrared HST observations by Kounkel et al. (2016). We therefore conclude that the two key properties

<caption>Table 2</caption>

Periodic Class I Sources Analyzed in this Study from Morales-Calderón et al. (2011)

<table><thead><tr><th>Source</th><th>Name</th><th><i>P</i> (day)</th><th>α (mag)</th><th>χ<sup>2</sup><sub>red</sub></th></tr></thead><tbody><tr><td>1</td><td>ISOY J053500.22-052409.2</td><td>7.29</td><td>0.39</td><td>3.12</td></tr><tr><td>2</td><td>ISOY J053516.33-052932.7</td><td>16.79</td><td>-0.19</td><td>2.72</td></tr><tr><td>3</td><td>ISOY J053520.72-051926.5</td><td>22.21</td><td>-0.06</td><td>4.09</td></tr><tr><td>4</td><td>ISOY J053522.10-051857.6</td><td>5.78</td><td>-0.08</td><td>2.27</td></tr><tr><td>5</td><td>ISOY J053524.73-051030.1</td><td>4.37</td><td>-0.08</td><td>1.91</td></tr><tr><td>6</td><td>ISOY J053525.18-050509.3</td><td>5.64</td><td>0.13</td><td>1.59</td></tr><tr><td>7</td><td>ISOY J053528.19-050341.2</td><td>23.52</td><td>-0.23</td><td>2.66</td></tr></tbody></table>

**Note.** Column (4) gives the amplitude $\alpha$ of the photometric variations associated with the spot, and Column (5) indicates the reduced chi squared value for the best-fit model (see Section 3.2).

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Color–color diagram of young stellar objects in Orion based on Spitzer/IRAC photometry. Gray plus signs and orange circles represent Class II sources and Class I sources, respectively, as identified by Megeath et al. (2012). Large black stars represent the sources considered in this study.</figcaption></figure>

defining our sample—periodicity and Class I nature—are truly associated with a single object in each source.

### 3.2. Fitting the Light Curves

We retrieved the photometric time series from Morales-Calderón et al. (2011) for the seven Class I periodic sources. Visual inspection revealed that, in at least several cases, the quasi-periodic variations that were initially identified through periodogram analysis occur within longer-timescale fluctuations. Morales-Calderón et al. (2011) suggest that these drifts reflect slow changes in the accretion rate or else self-shadowing in a warped disk. We remain agnostic about the nature of these long-timescale fluctuations but attempt to remove them in order to focus our analysis on the faster, manifestly periodic part of the signal. To this end, we subtract from the time series a linear trend obtained through least-squares fitting (see Figure 5). The amplitude of this linear variation over the 40 days of the Spitzer campaign ranges from 0.03 to 0.17 mag.

Once the photometric time series were detrended, we transformed them into phase-folded light curves using the periods proposed by Morales-Calderón et al. (2011). The scatter in the phased light curves is usually much larger than the quoted photometric precision (0.02–0.04 mag, typically), indicating that there is significant residual fluctuation from
<figure><img src="image_5.png" /><figcaption>Figure 4. SEDs of all seven objects in our sample (colored symbols), superimposed over a template SED derived from Taurus Class I sources (see Appendix A). Sources 5 and 7 have additional points, as they were included in the HOPS survey (Furlan et al. 2016).</figcaption></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>Figure 5. Original (black diamonds) and detrended (red triangles) light curve for Source 5. The black solid line is a linear fit to the original light curve, which is subtracted to obtain the detrended light curve. An arbitrary offset has been added to the detrended light curve for clarity.</figcaption></figure>

period to period even after detrending. Nonetheless, each light curve clearly departs from a flat line (see Figure 6), confirming the robustness of the original periodogram analysis.

We next assess whether the observed photometric variations can be ascribed to spots on the stellar surface. We focus on the 3.6 µm time series, but note that the same analysis of the 4.5 µm data leads to very similar results, albeit with an amplitude of the periodic signal that is slightly larger. As we will see, a simple model of a rotating star with a single spot on its surface does well to explain the observed, periodic variations. However, we acknowledge that the true source of variations could be multiple spots, or persistent disk features that somehow produce similar signals.

We assume that the spot covers a relatively small fraction of the stellar surface that does not overlap with the polar axis. We further assume that the spot itself can be idealized as an area of uniform brightness. The expected photometric signal introduced by the spot is then directly proportional to its projected area on the stellar surface. If the spot is small enough and disappears from the observer's view during each rotation period, this projected area is proportional to the cosine of the relative rotational phase $\phi-\phi_0$, where $\phi_0$ is the phase at which

the spot center crosses the stellar meridian. This idealized model ignores limb darkening, temperature gradients within spots, and the effects of large and/or multiple stellar spots that are commonly associated with young stellar objects (e.g., as seen in Doppler imaging; Skelly et al. 2008, 2009). Nonetheless, we believe that this model is sufficient for our limited goals—to assess whether the photometric variability is plausibly due to a rotating spot, and to provide a quantitative measure of the amplitude of that variability.

We thus proceed to perform a least-squares fit to each object's light curve, employing three parameters:

$$[3.6](\phi) = \begin{cases} [3.6]_0 + \alpha \cos [2\pi(\phi - \phi_0)] & \text{if } |\phi - \phi_0| \le 0.25 \\ [3.6]_0 & \text{otherwise.} \end{cases} \quad (4)$$

Here $[3.6]_0$ is the apparent 3.6 µm magnitude of the star when the spot is not in sight and $\alpha$ represents the amplitude of the spot-induced signal. In principle, the sign of $\alpha$ indicates the nature of the spot: a cold spot that makes the star temporarily fainter corresponds to a positive $\alpha$. Conversely, a hot spot induced by an accretion shock onto the stellar surface produces a negative $\alpha$.

We first explored parameter space through a fine sampling of $[3.6]_0$, $\alpha$, and $\phi_0$. When we considered only the quoted photometric errors in the data, even the best-fit model curves for each source had a relatively large $\chi^2$. The source of the problem was the residual, period-to-period fluctuations mentioned previously. To address this issue, we subtract the best-fit curve from the empirical data for each source and then calculate the “variability scatter” as the dispersion of the residuals. We quadratically add this scatter to the nominal photometric error to produce an “effective” uncertainty. This procedure assumes that the two errors are both random and uncorrelated with each other.

We then repeat our exploration of parameter space, but now using the enhanced, effective errors. Ideally, the resulting $\chi^2_{red}$ of each best-fit model should be close to unity. However, as seen in Tables 2 and 5, the actual $\chi^2_{red}$-values are still generally larger. We interpret this finding as an indication that the fluctuations are in fact correlated to some degree. In any case, we judge the model fits to be adequate for our purposes. Since
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Left: phased light curve for Source 6, an example of periodic dimming. The gray crosses and black diamonds represent the detrended and binned light curves, respectively, whereas the solid blue curve is the best-fitting “half-cosine” model. The binned light curve is shown for illustration only; the model is fit to the individual, detrended data points. The phase has been shifted so that the spot feature is centered around $\phi \approx 0.5$. Right: posterior probability distribution for the amplitude of the half-cosine model; the positive values of $\alpha$ indicate that the star gets temporarily fainter as the spot comes into view. Equivalent figures for the other sources are available in the electronic version of the journal.</figcaption></figure>

(The complete figure set (150 images) is available.)

$[3.6]_0$ and $\phi_0$ are not the focus of this study, we have marginalized over both these quantities to produce posterior probability distributions of $\alpha$ for each source. One representative light curve and $\alpha$-distribution is in Figure 6.

In all cases, the posterior distribution has one distinct peak. The corresponding $\alpha$-value of that peak is the amplitude of the best-fitting model; these values are listed in Table 2. The amplitudes range from 0.05 to 0.39 mag; weaker variations would probably remain undetectable given the *Spitzer* photometric precision and the intrinsic scatter around perfectly periodic signals displayed by our targets’ light curves. We note that in all cases models with $\alpha = 0$ are excluded by the probability distribution, confirming that there is a significant periodic signal in the light curve. Ostensibly, the prevalence of negative values of $\alpha$ for the best-fitting models suggests that the variability of most (but possibly not all) periodic Class I sources is due to accretion-generated hot spots.

### 3.3. Comparison with Class II and III Sources

It is instructive to compare the amplitudes of the periodically variable Class I sources in Orion with those of their Class II and III counterparts. For this purpose, we consider the sample of 150 variables whose periods were established by Morales-Calderón et al. (2011, see their Table 4). We assigned an evolutionary class to each object by locating it in the catalog of Megeath et al. (2012). If an object was not classified by Megeath et al. as a “protostar” (Class I) or “disk-bearing” (Class II), it was assumed to be Class III. Of the 150 variables, 7 are Class I, as previously noted, 65 are Class II, and 78 are Class III. Taking each object’s light curve from Morales-Calderón et al. (2011), we then performed our least-squares fit to establish the amplitude $\alpha$, in the same manner as for Class I sources. In almost all cases, the posterior distribution of $\alpha$-values shows one distinct peak. The distribution is more bimodal for a handful of sources with intrinsically weak variability. We list in Table B.1 the amplitude, period, and evolutionary class of all 143 objects.

While we were able to find least-squares fits to all Class II and III light curves, the limitations of our simple, half-cosine model became apparent with this larger sample. For instance, two of the Class II sources, ISOY J053605.95-050041.2 and ISOY J053518.03-052205.4, are eclipsing binaries (Morales-Calderón et al. 2011) and thus not well described by our model.

In addition, we would expect that Class III sources, corresponding to nonaccreting T Tauri stars, would have only cold spots ($\alpha > 0$). However, our fitting procedure finds that many have $\alpha < 0$. This apparent incongruity likely stems from the presence of multiple spots and departures from pure sinusoidal behavior even in the single-spot case. In summary, the *algebraic sign* of $\alpha$ cannot be interpreted literally. We maintain, however, that the *magnitude* of $\alpha$ still signifies quantitatively the degree of variability. Pursuing this line of reasoning, we find that the variability of Class I sources is much more similar to Class II objects than to those of Class III. Specifically, Class III sources generally have small values of $|\alpha|$, while Class I and II sources share a similar range in this quantity. The median values of $|\alpha|$ are 0.13, 0.11, 0.05 for Class I, II, and III, respectively. The left panel of Figure 7, which plots $\alpha$ against the period *P* for all three classes, shows the tendency. To confirm this impression more quantitatively, we performed a Kolmogorov–Smirnov test on the cumulative distributions of $|\alpha|$ for the three samples (see right panel of Figure 7). Table 3 lists the resulting values of $P_{KS}$ for the various pairs of samples. From these values, we can exclude the possibility that the Class I and III amplitudes were drawn from a single, underlying distribution. Similarly, the Class II and III sources do not share the same $|\alpha|$-distribution. On the other hand, periodic Class I and II sources have consistent amplitudes, indicating that the underlying physical objects could be identical. Previously, Wolk et al. (2015) found a similar result through observations of the GGD12-15 star-forming region (see their Figure 7). We ran a similar Kolmogorov–Smirnov test on the distribution of periods for each class of object. The results are shown in Table 3 and do
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Left: variability amplitude as a function of period for all Orion periodic sources. Classes I, II, and III are indicated by purple diamonds, blue triangles, and green squares, respectively. Right: cumulative histograms of the variability amplitude for each of the three types of sources. The solid black line denotes Class I, the dashed red line Class II, and the dotted-dashed blue line Class III.</figcaption></figure>

<table><caption>**Table 3**

P<sub>KS</sub> Values for Each Combination of Object as Obtained from a Kolmogorov-Smirnov Test on Their |α| and Period Distributions</caption><thead><tr><th>Object Categories</th><th>P<sub>KS</sub><br>|α|</th><th>P (day)</th></tr></thead><tbody><tr><td>Class I and Class II</td><td>0.33</td><td>0.21</td></tr><tr><td>Class I and Class III</td><td>0.003</td><td>0.39</td></tr><tr><td>Class II and Class III</td><td>3.1 10<sup>-5</sup></td><td>0.83</td></tr></tbody></table>

not indicate any significant difference in the distribution of periods between any two classes.

### 3.4. Envelope and Disk Models for Class I Light Curves

We have seen that both TSC and power-law ($\gamma = 1$) envelopes, in conjunction with circumstellar disks, can successfully reproduce the typical SED of a Class I source. We next examine whether endowing the central star with a rotating spot can yield a light curve at 3.6 µm that broadly matches the Orion observations. Following the current modeling of accreting, Class II sources, we consider the spot to be an optically thick, $10^4$ K layer of gas at the stellar surface (e.g., Hartmann et al. 2016). We place the spot at the star’s equator and select its surface area to be about 3% of the stellar surface so that the relative modulation at 3.6 µm from the star (with a disk but without any envelope) is 15%. This amplitude is motivated by the observations of periodic Class II variables in Orion (recall Figure 7). To generate our synthetic light curves for this asymmetric system, we use the 3D version of MCFOST. As in the previous, SED calculations, we set the observer’s line of sight at an inclination angle of 57°.

We examine first the light curve generated by a spotted star hosting a disk but without any envelope, i.e., a Class II source. In this case, thermal emission from the inner region of the disk overwhelms the stellar photosphere by a ratio of 4:1 at 3.6 µm. In addition, the area of the inner disk that directly faces the stellar hot spot is locally overheated, hence brighter. As a consequence, the light curve of such a system displays, in each period, two maxima in opposition of phase (see top panel in Figure 8). The first maximum, in phase with the stellar spot, has a much reduced amplitude compared to the diskless case. A

second maximum, due to local overheating of the disk by the hot spot, is observed in opposition of phase with the stellar spot, when this disk region is in the observer’s direct view.

We note parenthetically that the optical light curve of the star-disk system would exhibit a single maximum in phase with the spot, since disk emission is negligible in this regime. Cody et al. (2014) have observed this contrast between near-infrared and optical light curves in Class II sources, while Kesseli et al. (2016) have found a similar effect in simulations.

We now proceed to compute the light curve of a system containing a spotted star, a circumstellar disk, and an envelope, i.e., a full-fledged Class I system. Results from this analysis are shown in Figure 8 (middle and bottom panels).

Starting with the family of TSC models, we find that both the low-mass and best-fit models have only a very modest effect on the star-disk light curve, consistent with the fact that these envelopes are optically thin along the line of sight at 3.6 µm. The high-mass model, on the other hand, is sufficiently opaque to redistribute azimuthally the excess emission induced by the hot spot at the disk’s inner edge. As a result, the light curve only shows minor (<5%) variations that would only be marginally detectable through *Spitzer* monitoring. The $\gamma = 1$ power-law family of models is overall slightly more optically thick. The low-mass model is already sufficiently opaque to reduce the variability amplitude to ≈8%, still consistent with observations. On the other hand, the best-fit and high-mass models completely damp out any photometric variability.

It is well known that Class I sources drive polar outflows and jets that are believed to be responsible for the presence of polar cavities in their surrounding envelopes. This offers a different route for photons to escape the system. Besides their already-studied effects on the SED (e.g., Whitney et al. 2003; Furlan et al. 2016), one may wonder whether these cavities could alter significantly the light curve induced by a spotted central star. To test this possibility, we computed the light curves for additional models using our optimal TSC envelope structure but with a polar cavity carved out. The exact shapes of these cavities are not well established observationally, so we adopted two different morphologies, to explore the range of possibilities. Following Whitney et al. (2003), the first cavity shape we adopt is defined by the infall streamlines within the TSC model. The second cavity model has a simple conical shape. In both cases, we select

<figure><img src="image_10.png" /></figure>
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 8.** Model light curves for some models discussed in this study. The light curve for an isolated spotted star, indicated by black star symbols, is repeated in each panel. The light curve for the same star with a circumstellar disk is shown in the top left panel, whereas models of Class I sources using the TSC and power-law envelope models are displayed in the top right and bottom left panels, respectively. Finally, the bottom right panel illustrates the effect of introducing a polar cavity in the optimal TSC model.</figcaption></figure>

an outer opening angle of 20°, in line with scattered light images of Class I sources (Gramajo et al. 2010). Other curved cavity shapes have been explored in the past (Eisner et al. 2005; Furlan et al. 2008), but they are bracketed by our two options. In all cases, our chosen line of sight is at a higher inclination than the cavity opening. Thus, all the received light is reprocessed by the envelope. With too small a viewing angle, some light would travel directly from the star to the observer, and the resulting SED would be much more elevated at short wavelengths than in typical Class I sources.

The resulting light curves, presented in the bottom right panel of Figure 8, show that the effects of the cavity are modest, with the periodic signal maintaining a similar amplitude with and without a cavity. While photons can escape more easily from the system, most of them are not directed to an observer whose line of sight intercepts most of the envelope. Only the original signal, mildly attenuated by opacity through the envelope as we have shown before, reaches the observer. Finally, we note that the conical cavity reduces the amplitude of the signal more than the streamline cavity. In the latter case, the bounding streamline intersects the midplane 13 au from the central star. Since the region inside is evacuated, the optical depth along the observer's line of sight is diminished.

From this exploration, we conclude that some, but not all, of the SED-acceptable models of envelopes for Class I sources are compatible with the variations observed in the Orion Class I periodic sources. Specifically, the simultaneous constraints imposed by the objects' SED and their amplitude of variability imply that their envelopes must have $\tau_{3.6}$ of approximately 0.6–1.0 along the observer's line of sight. More opaque envelopes would reduce the photometric effects of the spotted star below their observed amplitude, while less opaque envelopes would result in too low a circumstellar extinction and thus too excessive near-infrared emission reaching the observer. For the TSC and $\gamma = 1$ power-law envelope models that we have considered in this study, the combined constraints from SEDs and light curves favor a total envelope mass of 0.025–0.05 $M_\odot$, i.e., an order of magnitude less than that of the central star.

## 4. Discussion

We now return to the fundamental question at the heart of our study, the physical nature of Class I sources. By studying the periodically variable ones, we first confirm that the central star in every case must be embedded in a three-dimensional circumstellar envelope. Unlike Class 0 sources, however, this structure must be of modest mass with respect to the star,
5%–10% of the stellar mass for our best models. Only then can the photometric variability, which certainly originates in the innermost region (given its timescale), be detected by an outside observer. An envelope that is too optically thick erases any photometric signal, both by reprocessing the emitted light to much longer wavelengths and by diffusively spreading any temperature anisotropy.

The central source that illuminates this moderately opaque envelope has photometric behavior that is consistent with Class II sources, which usually have disks, but not substantial envelopes. In other words, it appears that at least the periodically variable Class I sources in Orion are otherwise normal pre-main-sequence stars that happen to be embedded in moderately opaque envelopes. This picture corroborates previous findings that the spectroscopically determined effective temperatures, rotation rates, veiling, and line emission are similar in Class I and Class II sources (e.g., White & Hillenbrad 2004).

While mutually consistent, the dual constraints on the optical depth of Class I sources imposed by the SEDs and light curves ($0.5 \lesssim \tau_{3.6}^{\text{env}} \lesssim 4$ and $\tau_{3.6}^{\text{env}} \lesssim 1$, respectively) only partially overlap. Thus, the sample of Class I sources that display a periodic photometric signal is likely biased toward those with less opaque envelopes. Class I sources with nonperiodic light curves could be either (slightly) more embedded objects and/or aperiodic; the latter case is common among Class II sources (Morales-Calderón et al. 2011). Either way, we speculate, based on our analysis of the Orion periodic variables, that all Class I sources harbor normal pre-main sequence stars.

As for the internal structure of the envelope itself, our study provides only broad constraints. The widely used TSC envelope does yield a synthetic SED that roughly matches those of Class I sources. The fact has long been known (e.g., Whitney et al. 2003), but two caveats apply. First, we find that the envelope mass can be no greater than about $0.05 M_{\odot}$ to preserve the periodic signal in the 3.6 µm light curve. Within the framework of the underlying TSC model, this minimum mass corresponds to an envelope infall rate of $1.3 \times 10^{-6} M_{\odot} \text{yr}^{-1}$. However, we emphasize again that current observations inform us only about the envelope structure, and not its dynamics, and that this infall rate is in reality a proxy for the total envelope mass.

The geometric symmetry of Class I envelopes is still uncertain. In this regard, we have shown that inclusion of an outflow cavity has little impact on the observed SED or near-infrared light curve. Note that, even without an outflow cavity, TSC models are not spherically symmetric, since the density is lowered along the cloud's rotation axis. However, even this asymmetry is not essential for the envelopes to satisfy the two requirements of matching a typical Class I SED (see, e.g., Furlan et al. 2008, who modeled several Taurus Class I sources) and producing a satisfactory near-infrared light curve. We have illustrated this point by constructing much simpler envelope models that are spherically symmetric and have power-law density profiles, $\rho(r) \propto r^{-\gamma}$. Such models yield acceptable SEDs and light curves, as long as the envelope mass is relatively small and the power-law exponent $\gamma$ is roughly 1, the specific value we employed. Note that the azimuthally averaged density profile in the TSC model approaches, at small radii, a power law with $\gamma = 0.5$. Interestingly, the first imaging studies for Class 0 envelopes have found that pure power-law

envelope models yield a better fit than TSC models (see, e.g., Tobin et al. 2015; Persson et al. 2016).

More detailed study of Class I envelope structure will require spatially resolved mapping. Such observations need to be done at submillimeter or millimeter wavelengths, where the envelope is optically thin. In this regime, unfortunately, thermal emission from the cold outer regions of the circumstellar disk exceeds that from the envelope. Nonetheless, progress is being made. The PROSAC project used the SMA to disentangle envelope and disk masses in both Class 0 and I sources (Jorgensen et al. 2009).

If we are correct that the central stars in Class I sources are pre-main-sequence objects, then many will have hot spots on their surfaces, like their Class II counterparts. It is generally accepted that these are created by the transfer of mass from the inner regions of a circumstellar disk. The question then naturally arises as to whether a portion of the envelope is falling onto that disk to replenish it. Addressing this dynamical issue empirically will require spectral line observations of Class I envelopes. Such studies have also begun, both as part of the PROSAC project mentioned earlier and in targeted observations of individual sources (e.g., Leung et al. 2016). So far, the emphasis has been on the dynamics near outflow cavities; we hope studies will broaden to include the more general issue of cloud collapse.

Returning to our own work, it will be important for the future to increase the sample size of periodic Class I sources. The YSOVAR group has detected periodic variables in other star-forming regions (Günther et al. 2014; Poppenhaeger et al. 2015; Rebull et al. 2015; Wolk et al. 2015). It will be interesting to apply our analysis to these regions to expand the sample of Class I periodic sources.

To bolster the idea that the central objects are normal, pre-main-sequence objects, near-infrared spectroscopy would be helpful. Such observations would yield estimates of not only the effective temperature of the central objects but also the accretion rate on the stellar surface via continuum veiling and hydrogen line emission (e.g., Pa$\beta$, Br$\gamma$). Finally, Zeeman observations could reveal the magnetic field strength associated with these sources, to see if magnetic forces play a role in envelope structure.

The authors are grateful to Keivan Stassun for early discussions about this project, to Bill Herbst for input on the periodic variability of Class II and III sources, and to Christophe Pinte for his help in implementing and testing the addition of stellar spots in MCFOST. This research has made use of the NASA/IPAC Infrared Science Archive, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation.

## Appendix A
### A Class I Template SED

It is necessary to assemble the observed SED of a Class I source in order to evaluate which of our model SEDs are representative of such objects. The SED of the targets studied in this work only extends out to 24 µm (and often only 8 µm). While this is sufficient to assess their embedded nature through
<figure><img src="image_15.png" /><figcaption>**Figure 9.** Left: luminosity-normalized SEDs of Taurus Class I sources. Gray crosses represent individual photometric measurements, whereas the red diamonds and error bars indicate the median and 16th-to-84th percentile range for wavelengths with observations of at least 10 objects (out of the sample of 21). Right: luminosity-normalized SEDs of Orion Class I sources (blue triangles and error bars indicate the median and 16th-to-84th percentile range, respectively). Red diamonds and error bars reproduce the Taurus Class I template.</figcaption></figure>

their near- to mid-infrared spectral index, the model SEDs must be compared over a much broader range of wavelengths, spanning from the optical/near-infrared to the millimeter regime. Furthermore, there is a rich diversity of Class I SEDs, possibly due to viewing angle effects and an intrinsic diversity in envelope size and mass. To draw a representative ("template") SED, it is therefore necessary to consider a larger sample of Class I sources and to use their average and dispersion. Similar efforts have been conducted for Class II sources in the past (e.g., Mathews et al. 2013).

To ensure the richest and most uniform data set possible, we have focused on the population of Class I sources in the nearby Taurus-Auriga star-forming region, arguably the best-studied young stellar population. We based our sample on the survey conducted by Furlan et al. (2008). From this initial sample, we excluded CoKu Tau 1, DG Tau B, HH 30, IRAS 04016+2610, IRAS 04158+2805, IRAS 04248+2612, and IRAS 04302+2247, in which a nearly edge-on circumstellar disk strongly affects the system's SED (Burrows et al. 1996; Padgett et al. 1999; Glauser et al. 2008). Indeed, in some extreme cases (HH 30, IRAS 04158+2805), an edge-on disk can make a Class II source erroneously appear like a Class I embedded object if one only focuses on the near- to mid-infrared spectral index.

For the 21 bona fide Class I sources in the sample, we adopted the broadband photometry presented in Furlan et al. (2008) and complemented it with photometric measurements from the *AKARI* and *WISE* all-sky surveys whenever available. We also included *Herschel* photometry for a few individual sources from Howard et al. (2013), Bulger et al. (2014), and Ribas et al. (2017) in cases where far-infrared observations were lacking. The resulting photometric database contains 455 individual flux density measurements for 21 objects at wavelengths ranging from 0.35 µm to 1.3 mm.

The range of bolometric luminosities for this sample is very broad, from 0.25 to 28 $L_{\odot}$ (Furlan et al. 2008), such that distribution of the objects' flux densities at a given wavelength does not just reflect the diversity in the shape of Class I SEDs but also the range of their luminosities. To sidestep this problem, we transformed each object's SED to the unitless

quantity $\lambda F_{\lambda} 4\pi d^2/L_{\text{bol}}$ (where $d = 140$ pc), which effectively indicates the proportion of the object light output as a function of wavelength, i.e., the global shape of the SED.

We then grouped all luminosity-normalized measurements by observing wavelengths. At each wavelength for which more than 10 objects have a measured flux density (most notably, this limits our analysis to $\lambda \ge 1$ µm), we then compute the median and the 16th and 84th percentiles. These estimates define our “template (Taurus) Class I SED” and associated uncertainties. Note that the latter rather represent the range of possible SED shapes, or intrinsic scatter within the population, rather than the actual measurement uncertainties. The template SED, along with all individual measurements, is shown in the left panel of Figure 9 and listed in Table 4. The dispersion between objects is markedly larger at the longest and, especially, the shortest wavelengths, whereas most objects are within a factor of ≈2 of the median SED in the mid- to far-infrared regime. This is a natural consequence of our normalization scheme and of the fact that this is the regime where most of the luminosity of a Class I is radiated away.

While the Taurus Class I sources are the best-studied sample available to build such a template, the sample is relatively small and may not be representative of Class I sources in other star-forming regions. To address this issue, we followed the same procedure for the sample of 219 Orion Class I and flat-spectrum sources from Furlan et al. (2016). As illustrated in the right panel of Figure 9, the general structure and ranges of the resulting Orion Class I template are consistent with the Taurus template. The only marginal difference between the two templates is in the submillimeter regime, where the Orion template is higher than the Taurus one by a factor of a few. Most likely, this is a consequence of the enhanced crowding of the Orion region relative to Taurus, which is further compounded by confusion with cloud emission. These effects lead to a Malmquist bias whereby only the brightest sources can be detected in Orion. Indeed, in the submillimeter, only 20%–30% of all Orion sources are detected and included in the template, whereas in Taurus the corresponding fraction is roughly 70%. Given the similarity of the two templates, we

<figure><img src="image_16.png" /></figure>
<table><caption>**Table 4**
Luminosity-normalized Template SED for Class I Sources in Taurus</caption><thead><tr><th rowspan="2">λ</th><th colspan="3">λF<sub>λ</sub> 4πd<sup>2</sup>/L<sub>bol</sub></th><th rowspan="2">N<sub>obj</sub></th><th rowspan="2">Filter/Mission</th></tr><tr><th>16th Percentile</th><th>Median</th><th>84th Percentile</th></tr></thead><tbody><tr><td>1.25</td><td>0.000281</td><td>0.00291</td><td>0.0308</td><td>19</td><td>2MASS J</td></tr><tr><td>1.65</td><td>0.00188</td><td>0.0355</td><td>0.0846</td><td>20</td><td>2MASS H</td></tr><tr><td>2.20</td><td>0.00763</td><td>0.0536</td><td>0.160</td><td>21</td><td>2MASS K</td></tr><tr><td>3.6</td><td>0.0192</td><td>0.0861</td><td>0.318</td><td>19</td><td>Spitzer/IRAC</td></tr><tr><td>4.65</td><td>0.0182</td><td>0.108</td><td>0.219</td><td>17</td><td>Spitzer/IRAC</td></tr><tr><td>5.8</td><td>0.0105</td><td>0.137</td><td>0.267</td><td>13</td><td>Spitzer/IRAC</td></tr><tr><td>8.</td><td>0.00593</td><td>0.0915</td><td>0.210</td><td>11</td><td>Spitzer/IRAC</td></tr><tr><td>9.</td><td>0.104</td><td>0.143</td><td>0.317</td><td>16</td><td>AKARI</td></tr><tr><td>10.5</td><td>0.0430</td><td>0.118</td><td>0.193</td><td>15</td><td>N</td></tr><tr><td>12.</td><td>0.0733</td><td>0.190</td><td>0.281</td><td>20</td><td>IRAS</td></tr><tr><td>18.</td><td>0.148</td><td>0.225</td><td>0.501</td><td>17</td><td>AKARI</td></tr><tr><td>22.1</td><td>0.159</td><td>0.249</td><td>0.466</td><td>18</td><td>WISE</td></tr><tr><td>25.</td><td>0.130</td><td>0.295</td><td>0.447</td><td>18</td><td>IRAS</td></tr><tr><td>60.</td><td>0.243</td><td>0.373</td><td>0.437</td><td>18</td><td>IRAS</td></tr><tr><td>65.</td><td>0.179</td><td>0.293</td><td>0.400</td><td>14</td><td>AKARI</td></tr><tr><td>90.</td><td>0.120</td><td>0.242</td><td>0.432</td><td>15</td><td>AKARI</td></tr><tr><td>100.</td><td>0.124</td><td>0.320</td><td>0.456</td><td>14</td><td>IRAS</td></tr><tr><td>140.</td><td>0.0874</td><td>0.168</td><td>0.394</td><td>15</td><td>AKARI</td></tr><tr><td>160.</td><td>0.0801</td><td>0.157</td><td>0.491</td><td>14</td><td>AKARI</td></tr><tr><td>350.</td><td>0.00705</td><td>0.0312</td><td>0.0488</td><td>10</td><td></td></tr><tr><td>450.</td><td>0.00219</td><td>0.00736</td><td>0.0249</td><td>17</td><td></td></tr><tr><td>800.</td><td>0.000380</td><td>0.000903</td><td>0.00218</td><td>14</td><td></td></tr><tr><td>850.</td><td>0.000347</td><td>0.000914</td><td>0.00207</td><td>17</td><td></td></tr><tr><td>1100.</td><td>0.000136</td><td>0.000275</td><td>0.000589</td><td>13</td><td></td></tr><tr><td>1300.</td><td>0.000121</td><td>0.000273</td><td>0.00127</td><td>17</td><td></td></tr></tbody></table>

**Note.** The fifth column indicates the number of objects that contribute to the estimate of the median and the 16th and 84th percentiles of the distribution at each wavelength. Only wavelengths with $N_{obj} \ge 10$ are listed.

adopt for our study the more uniform and complete Taurus template.

## Appendix B
150 Variable Sources in Orion

We considered the 150 variable objects in Orion that had newly identified periods in Morales-Calderón et al. (2011, see their Table 4) to place the amplitude of variability of the Class I sources studied here in their broader context. Roughly half of these 150 objects were assigned the classification given to them in Megeath et al. (2012), where we considered “protostars” to be Class I (these are the same targets as studied in Section 3) and “disk-bearing” objects to be Class II. We then classified the remaining objects that were not assessed in Megeath et al. (2012) as Class III. Following the procedure described in Section 3.2, the light curves from these 150 sources were fit to a single spotted model using the method of least squares. The

<table><caption>**Table 5**
Amplitude of Variation Pulled from the Best Fit to the Object's Light Curve, Listed with the Object's Name, Class, and the Reduced Chi Square Value from the Best-fitting Light Curve</caption><thead><tr><th>Object Name</th><th>Class</th><th>P (day)</th><th>α</th><th>χ<sup>2</sup><sub>red</sub></th></tr></thead><tbody><tr><td>ISOY J053427.70-053155.4</td><td>III</td><td>3.37</td><td>-0.05</td><td>2.63</td></tr><tr><td>ISOY J053434.88-044243.5</td><td>III</td><td>9.36</td><td>0.05</td><td>3.03</td></tr><tr><td>ISOY J053435.15-053210.4</td><td>III</td><td>5.18</td><td>-0.01</td><td>1.14</td></tr><tr><td>ISOY J053435.98-045218.0</td><td>II</td><td>2.34</td><td>0.11</td><td>2.67</td></tr><tr><td>ISOY J053439.11-053839.2</td><td>III</td><td>2.41</td><td>0.05</td><td>2.29</td></tr><tr><td>ISOY J053439.76-052425.6</td><td>III</td><td>16.97</td><td>0.05</td><td>2.65</td></tr><tr><td>ISOY J053441.18-054611.7</td><td>II</td><td>8.02</td><td>0.11</td><td>1.98</td></tr><tr><td>ISOY J053441.71-052653.0</td><td>II</td><td>8.47</td><td>0.05</td><td>1.73</td></tr><tr><td>ISOY J053441.87-055502.6</td><td>II</td><td>1.66</td><td>-0.06</td><td>2.55</td></tr><tr><td>ISOY J053442.74-052837.6</td><td>II</td><td>6.15</td><td>-0.12</td><td>3.26</td></tr><tr><td>ISOY J053445.00-045559.2</td><td>II</td><td>3.40</td><td>-0.07</td><td>1.71</td></tr><tr><td>ISOY J053447.64-054350.7</td><td>III</td><td>4.27</td><td>0.05</td><td>2.00</td></tr><tr><td>ISOY J053447.74-052632.1</td><td>II</td><td>1.84</td><td>0.11</td><td>1.99</td></tr><tr><td>ISOY J053447.85-044228.9</td><td>III</td><td>8.36</td><td>-0.08</td><td>2.08</td></tr><tr><td>ISOY J053449.57-052903.4</td><td>II</td><td>14.73</td><td>-0.13</td><td>1.95</td></tr><tr><td>ISOY J053450.66-050407.7</td><td>III</td><td>9.24</td><td>0.02</td><td>1.06</td></tr><tr><td>ISOY J053450.72-045836.8</td><td>II</td><td>6.32</td><td>0.12</td><td>3.94</td></tr><tr><td>ISOY J053452.22-045102.7</td><td>III</td><td>5.15</td><td>-0.05</td><td>2.16</td></tr><tr><td>ISOY J053452.49-044940.4</td><td>III</td><td>6.25</td><td>-0.01</td><td>1.23</td></tr><tr><td>ISOY J053452.60-051536.6</td><td>III</td><td>2.35</td><td>0.13</td><td>5.31</td></tr></tbody></table>

**Note.** Objects marked with † are eclipsing binary candidates. Only the first 20 entries in the table are shown here; the full table is available electronically.
(This table is available in its entirety in machine-readable form.)

results from our light-curve analysis, specifically the variability amplitude α, are presented in Table 5.

## ORCID iDs

Gaspard Duchene https://orcid.org/0000-0002-5092-6464

## References

Allard, F., Homeier, D., & Freytag, B. 2012, RSPTA, 370, 2765

André, P., Ward-Thompson, D., & Barsony, M. 1993, ApJ, 406, 122

Beck, T. L. 2007, AJ, 133, 1673

Bulger, J., Patience, J., Ward-Duong, K., et al. 2014, A&A, 570, A29

Burrows, C. J., Stapelfeldt, K. R., Watson, A. M., et al. 1996, ApJ, 473, 437

Cassen, P., & Moosman, A. 1981, Icar, 48, 353

Chiang, H.-F., Looney, L.-W., & Tobin, J. J. 2012, ApJ, 756, 168

Cody, A. M., Stauffer, J., Baglin, A., et al. 2014, AJ, 147, 82

Connelley, M. S., & Greene, T. P. 2010, AJ, 140, 1214

Covey, K. R., Greene, T. P., Doppmann, G. W., & Lada, C. J. 2005, AJ, 129, 2765

Doppmann, G. W., Greene, T. P., Covey, K. R., et al. 2005, AJ, 130, 1145

Eisner, J. A., Hillenbrand, L. A., Carpenter, J. M., et al. 2005, ApJ, 635, 396

Evans, N. J., II, Dunham, M. M., Jørgenson, J. K., et al. 2009, ApJ, 181, 321

Flaherty, K. M., Pipher, J. L., Megeath, S. T., et al. 2007, ApJ, 663, 1069

Furlan, E., Fischer, W. J., Ali, B., et al. 2016, ApJS, 224, 5

Furlan, E., McClure, M., Calvet, N., et al. 2008, ApJS, 176, 184

Glauser, A. M., Ménard, F., Pinte, C., et al. 2008, A&A, 485, 531

Gramajo, L. V., Whitney, B. A., Gómez, M., et al. 2010, AJ, 139, 2504
Günther, H. M., Cody, A. M., Covey, K. R., et al. 2014, AJ, 148, 122

Haisch, K. E., Lada, E. A., & Lada, C. J. 2001, ApJ, 553, 153

Hartmann, L., Herczeg, G., & Calvet, N. 2016, ARA&A, 54, 135

Howard, C. D., Sandell, G., Vacca, W. D., et al. 2013, ApJ, 776, 21

Jorgensen, J. K., van Dishoeck, E. F., Visser, R., et al. 2009, A&A, 507, 861

Kenyon, S. J., & Hartmann, L. 1995, ApJS, 101, 117

Kesseli, A. Y., Petkova, M. A., Wood, K., et al. 2016, ApJ, 828, 42

Koresko, C. D., Herbst, T. M., & Leinert, C. 1997, ApJ, 480, 741

Kounkel, M., Megeath, S. T., Poteet, C. A., Fischer, W. J., & Hartmann, L. 2016, ApJ, 821, 52

Leung, G. Y. C., Lim, J., & Takakuwa, S. 2016, ApJ, 833, 55

Mathews, G. S., Pinte, C., Duchêne, G., Williams, J. P., & Ménard, F. 2013, A&A, 558, A66

Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425

Megeath, S. T., Gutermuth, R., Muzerolle, J., et al. 2012, AJ, 144, 192

Morales-Calderón, M., Stauffer, J. R., Hillenbrand, L. A., et al. 2011, ApJ, 733, 50

Muzerolle, J., Hartman, L., & Calvet, N. 1998, AJ, 116, 2965

Myers, P. C., & Ladd, E. F. 1993, ApJ, 413, 47

Padgett, D. L., Brandner, W., Stapelfeldt, K. R., et al. 1999, AJ, 117, 1490

Persson, M. V., Harsono, D., Tobin, J. J., et al. 2016, A&A, 590, A33

Pinte, C., Harries, T. J., Min, M., et al. 2009, A&A, 498, 967

Pinte, C., Ménard, F., Duchêne, G., & Bastien, P. 2006, A&A, 459, 797

Poppenhaeger, K., Cody, A. M., Covey, K. R., et al. 2015, AJ, 150, 118

Prato, L., Lockhart, K. E., Johns-Krull, C. M., et al. 2009, AJ, 137, 3931

Rebull, L. M., Stauffer, J. R., Cody, A. M., et al. 2015, AJ, 150, 175

Ribas, Á., Espaillat, C. C., Macías, E., et al. 2017, ApJ, 849, 63

Robitaille, T. P., Whitney, B. A., Indebetouw, R., Wood, K., & Denzmore, P. 2006, ApJS, 167, 256

Shu, F. H. 1977, ApJ, 214, 488

Skelly, M. B., Unruh, Y. C., Barnes, J. R., et al. 2009, MNRAS, 399, 1829

Skelly, M. B., Unruh, Y. C., Collier Cameron, A., et al. 2008, MNRAS, 385, 708

Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163

Tobin, J. J., Stutz, A. M., Megeath, S. T., et al. 2015, ApJ, 798, 128

White, R. J., & Hillenbrad, L. A. 2004, ApJ, 616, 998

Whitney, B. A., Wood, K., Bjorkman, J. E., & Wolff, M. J. 2003, ApJ, 591, 1049

Wolk, S. J., Günther, H. M., Poppenhaeger, K., et al. 2015, AJ, 150, 145

Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868