The JCMT Transient Survey: Four-year Summary of Monitoring the Submillimeter Variability of Protostars

We first examine secular variability, i.e., changes that take place across the entire time coverage. We determine the amplitudes and timescales associated with the observed brightness fluctuations, and then use these measurements to classify the types of variability. Furthermore, these estimates for the timescales involved in the observed variations provide a direct link to the underlying physical process responsible for the variations (e.g., viscously evolving disks, Keplerian orbit times).

We apply an LSP (Lomb 1976; Scargle 1989; VanderPlas 2018) analysis in order to uncover and quantify the relevant timescales and amplitudes of secular variability observed in the submillimeter. The LSP measures the likelihood that an observed light curve is fit well by a sinusoidal function with a set of parameters: mean brightness, amplitude, frequency, and phase. The statistical power measures the quality of the best sinusoidal fit as a function of frequency and builds the periodogram of the observed light curve. We perform the LSP analysis using the LombScargle task in the timeseries package of astropy (Astropy Collaboration et al. 2013) and analytically modify the output to fit our purpose. The details of the LSP output and our modifications are described below.

In the LSP analysis, the statistical power P_A of the "A" hypothesis is defined as the fractional improvement over the non-varying (or static) "N" hypothesis:

$$\begin{eqnarray}&&{P}_{{\rm{A}}}=\displaystyle \frac{{\chi }_{N}^{2}-{\chi }_{{\rm{A}}}^{2}}{{\chi }_{N}^{2}},\end{eqnarray} \tag{ 1 }$$

where ${\chi }_{X}^{2}$ denotes the chi-squared value of the subscripted hypothesis. Thus, the value of P_A rises as the goodness of fit of the sinusoidal function improves on the non-varying hypothesis; that is, ${\chi }_{{\rm{A}}}^{2}$ reduces with respect to ${\chi }_{N}^{2}$. Following standard practice, we consider the frequency at which P_A reaches its maximum value to denote the best fit.

We next determine the False Alarm Probability (FAP) of this best fit to validate the detection. The FAP measures the probability that the sinusoidal fit is not real. The FAP is calculated using the best-fit χ² (Baluev 2008). That is, the FAP of the best-fit "A" hypothesis, FAP_A, is defined as

$$\begin{eqnarray}&&{{FAP}}_{{\rm{A}}}={\left(1-{P}_{{\rm{A}}}\right)}^{{N}_{{\rm{f}}}/2}.\end{eqnarray} \tag{ 2 }$$

The FAP decreases with increasing degrees of freedom in the fit N_f = (N_epoch − D_A), where N_f is determined from the number of epochs observed, N_epoch, minus the number of parameters in the model, D_A. The FAP_A also reduces as the fractional improvement of the sinusoidal fit P_A increases.

Following the traditional LSP approach, the FAP of the best hypothesis must also take into account the total number of hypotheses tested. Thus, we define FAP_LSP as

$$\begin{eqnarray}&&{{FAP}}_{\mathrm{LSP}}=1-{\left(1-{{FAP}}_{{\rm{A}}}\right)}^{{N}_{\mathrm{freq}}},\end{eqnarray} \tag{ 3 }$$

where N_freq is the total number of frequencies tested. Furthermore, since we are only interested in cases where FAP_A ≪ 1, this can be significantly simplified to

$$\begin{eqnarray}&&{{FAP}}_{\mathrm{LSP}}\cong {{FAP}}_{{\rm{A}}}\times {N}_{\mathrm{freq}}.\end{eqnarray} \tag{ 4 }$$

Most protostellar sources are best fit by low-frequency, long-timescale, sinusoids. Very short-period sinusoids always result in poor fits to these protostellar light curves. The tested frequency range and step is calculated using the time baseline and the number of observed epochs. Using our 4 yr time baseline and 35 epochs, the tested frequencies are between 0.025 yr⁻¹ (once in 40 yr) and 22 yr⁻¹ (almost once per two weeks) with 0.05 yr⁻¹ steps. The traditional FAP_LSP calculated over the full tested frequency range, which includes very high frequencies, will provide a systematic overestimate of the false alarm probability. We therefore modify the FAP_LSP to consider only periods longer than the best-fit period, with ${{FAP}}_{\mathrm{Mod}}$ given by

$$\begin{eqnarray}&&{{FAP}}_{\mathrm{Mod}}\cong {{FAP}}_{{\rm{A}}}\times {N}_{\lt \mathrm{freq}},\end{eqnarray} \tag{ 5 }$$

where N_<freq is the number of tested periods that are longer than or equal to the best-fit period (lower than the best-fit frequency). Thus, ${{FAP}}_{\mathrm{Mod}}$ becomes equivalent to FAP_LSP for sources with best-fit periods similar to the highest frequency tested.

Given that our LSP analysis uncovers primarily long-period sinusoids, often with periods greater than the 4 yr observing window, we also apply a linear least-squares fitting method to quantify the best-fit linear solution. The statistical power, P_Lin, and the FAP of the linear best fit, FAP_Lin, are calculated following Equations (1) and (2).

We demonstrate our dual secular analysis procedure in Figure 1, where the bottom panel plots the 4 yr light curve for a known variable (EC 53 in Serpens Main; see also Lee et al. 2020b). We first perform the LSP analysis and reproduce the periodogram in the upper left panel. We next diagnose the various FAPs derived at the peak of the periodogram, denoted by the vertical orange line. For this source, the periodogram peak is found at a frequency of 0.0017 day⁻¹, equivalent to a 1.61 yr period, the single-frequency false alarm probability is less than 10⁻¹¹, and both the LSP and modified FAPs are less than 10⁻⁸, denoting that the sinusoidal fitting is highly faithful. Finally, we also determine the best-fit linear form (dashed orange line in upper panel) and calculate a large linear FAP, greater than 0.5, denoting the linear fit has a 50% likelihood of being a false alarm.

Zoom In Zoom Out Reset image size

In Figure 2, we show examples of sources with light curves fit robustly by either periodic or linear secular functions. The top panel shows a periodic light curve with short period compared to our time coverage (4 yr). The second panel from the top shows a light curve robustly fit by a sinusoid with a period comparable to 4 yr, but not by a linear function. The bottom two panels show light curves fit well by both sinusoidal and linear functions, where the sinusoidal period is much longer than the time baseline. We specify the classification of the found variables and analyze their properties in Section 3.3.

Zoom In Zoom Out Reset image size

To confirm the robustness of our method, in Figure 3 we compare the FAPs of the two fitting functions, sinusoidal and linear, for each of the 295 submillimeter sources in our sample. In the figure, we color-code the sources by type: protostar (red), disk (blue), and starless (gray). The starless cores and disk sources occupy the lower left corner of the plot, where there is no robust evidence for secular variability (i.e., FAPs > 0.1%).

Zoom In Zoom Out Reset image size

Eighteen protostellar sources are found significantly outside the lower left corner of Figure 3. Furthermore, all 14 protostellar sources that have FAP_Lin < 0.001 also have ${{FAP}}_{\mathrm{Mod}}\lt 0.001$ and lie along a diagonal line in the figure. This is expected, given that long-period sinusoids are almost indistinguishable from linear functions over limited timescales. In the figure, however, there are an additional four protostellar sources with ${{FAP}}_{\mathrm{Mod}}\lt 0.001$ for which a linear fit is poor due to the observed periodicity or curvature of the light curve. While EC 53 (see Figure 1) shows clear periodic variability, the other three sources, HOPS 315, HOPS 373, and SMM 10, reveal clear curved trends with little overall slope and have best-fit periods of < 8 yr; see the light curves of HOPS 315 (Figure B1.9), HOPS 373 (Figure 2), and SMM 10 (Figure B1.5). We discuss further all of the revealed secular variables in Section 3.3.

It is important to note that the LSP determination of a robust sinusoid fit does not necessarily imply that the underlying source variability is truly periodic, especially for the majority of our uncovered sources (16 out of 18) requiring long periods, on the order of or greater than the 4 yr observing time window. Nevertheless, the best-fit parameters provide a useful, quantifiable estimate of the underlying timescale and amplitude for observed light-curve variations. We stress, however, that for periods longer than 4 years, these are extrapolated estimates of the period and amplitude of the observed light-curve event and not an indication of underlying periodicity. Additionally, for the longest robust sinusoidal fits, the periods of the amplitudes are extremely uncertain and the linear slopes are more appropriate quantifiable measures of the variability. Thus, in general we will use sinusoidal fits only for sources with derived periods less than ∼15 years.

Along with the long-term secular variability, sources may vary irregularly, with timescales of a few months or shorter. Two types of stochasticity are anticipated: (1) long-term chaotic brightness variations resulting in an increase in the epoch-to-epoch brightness variability without an appreciable linear or periodic trend, perhaps due to short-timescale instabilities near the disk/star interface; (2) rare, individual epoch brightening or dimming events that robustly stand out above the noise. We explore both of these stochastic variability possibilities across all observed epochs for each of the 295 submillimeter sources in our sample.

We first consider the rare, individual epoch events. The expected uncertainty of each peak brightness measurement of a given submillimeter source at 850 μm is defined as the fiducial uncertainty, σ_fid, and calculated as

$$\begin{eqnarray}&&{\sigma }_{\mathrm{fid}}=\sqrt{{0.014}^{2}+{\left(0.02\times \bar{F}\right)}^{2}}\ \mathrm{Jy}\,{\mathrm{beam}}^{-1},\end{eqnarray} \tag{ 6 }$$

where $\bar{F}$ is the mean peak brightness of the source over all observations. The other two terms in Equation (6) are the relative flux calibration error for a given map epoch, 2%, and the typical 850 μm rms noise, 0.014 Jy beam⁻¹. These values have been empirically measured for the JCMT Transient Survey by Mairs et al. (2017a), and the fiducial uncertainty is described in detail by Johnstone et al. (2018). To mitigate the effect of outlier measurements on the mean peak brightness, we determine the mean peak brightness after excluding the brightest and faintest 10% of measurements from each source, ${\bar{F}}_{80 \% }$, and we use this value in place of $\bar{F}$ in Equation (6). Furthermore, for the 18 secular variables uncovered in Section 3.1, we subtract the best-fit sinusoidal function from their light curves before determining ${\bar{F}}_{80 \% }$.

We first analyze the light curves of all 295 submillimeter sources in our sample for individual outlier events. For each measurement of each source, including the brightest and faintest 10%, we determine the residual peak brightness after subtracting ${\bar{F}}_{80 \% }$. We then bring these residual measurements to a standard scale by dividing by the expected measurement uncertainty σ_fid. Figure 4 presents a histogram over all these σ_fid scaled residual fluxes (light), and after excluding the 18 secular variables (dark).

Zoom In Zoom Out Reset image size

The shape of the distribution follows well the expected normal distribution with an standard deviation fit close to unity, σ = 1.10, indicating that our uncertainty estimates are appropriate. The residuals with absolute values greater than 3, however, are somewhat overpopulated, perhaps indicating slight non-Gaussian shape in the noise characteristics. In addition, there are six highly stochastic events with absolute residuals greater than 5, indicated by the red vertical lines in the figure. These six candidate stochastic events (Table 2) are measured across six different sources, five of which are known protostars,

Table 2. Candidate Stochastic Events

Region	ID	Known Name	${\bar{{\rm{F}}}}_{80 \% }$($\bar{{\rm{F}}}$) a	${\sigma }_{\max }$	Fmax/${\bar{{\rm{F}}}}_{80 \% }$	Date	Robust?
				(Jy beam−1)	(σfid)		(light curve)
NGC 2068	J054631.0−000232	HOPS 373	1.25 (1.26)	5.34 b	1.12	2019 Sep 26	Yes (Figure 2)
OMC2/3	J053522.4−050111	HOPS 88	2.50 (2.50)	6.80	1.14	2016 Feb 29	No/Focus
Serpens south	J183002.6−020248	CARMA 3	2.19 (2.20)	6.43	1.13	2019 Nov 2	No/Focus
OMC2/3	J053527.4050929−	HOPS 370	2.62 (2.64)	5.79	1.11	2017 Apr 21	No/Focus
NGC 2068	J054645.6 + 000719	LBS 11	0.32 (0.32)	5.48	1.26	2016 Apr 27	No/Extended
OMC2/3	J053524.2−050932c	starless	0.34 (0.34)	−5.01		2018 Nov 2	No/Extended

Notes.

^a The values in the bracket denote the mean peak brightness using the whole data points (see Section 3.2.) ^b Value obtained after subtracting the best-fit sinusoidal function.

Download table as:  ASCIITypeset image

Two of the six candidate stochastic events are associated with faint submillimeter sources, ${\bar{F}}_{80 \% }\sim 0.3$ Jy beam⁻¹, for which the fiducial uncertainty is dominated by measurement error (Equation (6)). This sample includes the only non-protostellar potential variable source in our sample, and also the only source to apparently dim for a single epoch. Careful examination of the residual images, after subtracting the co-add over all epochs from the epoch in which the event occurred, reveals that for both these candidates the residual emission is extended significantly beyond that of a point source. These results strongly suggest that for these two epochs the image reconstruction procedure has introduced spurious large-scale features, not unexpected in submillimeter map reconstructions (for details, see Mairs et al. 2015, 2017a). We therefore drop these two candidates as potential stochastic events.

Three of our candidate stochastic events are bright, ${\bar{F}}_{80 \% }\,\sim 2$ Jy beam⁻¹, but reside nearby, r < 30'', vastly brighter sources, ${\bar{F}}_{80 \% }\gt 5$ Jy beam⁻¹. For each of these candidates, the residual images for the epochs of interest show clear structure in the residual beam pattern of all the brighter sources within the map, indicating slight focus issues during the observation. The focus issue is seen to clearly produce excess emission within r ∼ 30'' around the brighter sources, and thus the measurement uncertainties for the three candidate stochastic events are significantly underestimated at these particular epochs. We therefore drop these three candidates as potential stochastic events.

Finally, one of our protostellar sources, HOPS 373 (see Figure 2 and Table 2), is detected in both the secular and stochastic variability analyses. We therefore note that the stochastic event detected in HOPS 373 appears to be part of a longer-timescale burst-like event. Comparison of the residual map for the epoch in which the >5σ stochastic event occurred reveals a strong point-like significant peak at the location of HOPS 373, assuring the robustness of the detection. Given that this brightening event occurs over more than one epoch, ³⁶ we categorize this source as a secular variable in the rest of this work.

Mairs et al. (2019) performed a separate single-epoch transient source analysis for those areas within our monitored star-forming regions where there are no bright submillimeter sources. That investigation uncovered a prominent stochastic submillimeter variable, coincident with JW 566 in OMC2/3. JW 566 is a T Tauri star and shows a brightening by 500 mJy in a single epoch and a ∼50% decline in its peak brightness during the half-hour observation. The mean brightness of JW 566 over all epochs lies below the threshold used for this paper, however, and therefore JW 566 is not identified by, or included in, our analysis.

Having found no robust rare, individual epoch events within our sample of bright sources, we next consider the measured peak brightness uncertainty over the full light curve, normalized by the expected fiducial uncertainty, as a determination of the long-term stochastic nature of each of our 295 submillimeter sources. Here, we are searching for sources with a significantly larger spread in brightness measurements compared with the known noise properties of the measurements. Figure 5 plots this dispersion against mean peak brightness, with color again indicating source type. The y-axis plots ${\chi }_{80 \% }^{2}$, the chi-square obtained by using only the flux measurements that went into calculating ${\bar{F}}_{80 \% }$. Thus, any rare significant outlier events are not included in this analysis. The eighteen known secular variables are marked with black circles, and HOPS 373 with the individual stochastic event is marked with cross. The left panel presents the calculation before removing the best-fit sinusoid from the secular variables, whereas the right panel shows the result after subtracting the secular trend.

Zoom In Zoom Out Reset image size

Two results are immediately evident in Figure 5. First, all sources in the left panel with anomalously large variations in the measured peak brightness over their full light curves, ${\chi }_{80 \% }^{2}/{N}_{80 \% }\gt 1.8$ (corresponding to an expected χ²/N ∼ 4 when measured over a full Gaussian distribution), are secular variables for which the large deviations are significantly removed by subtracting the secular fit (right panel). Thus, given the present survey sensitivity, there are no sources with ongoing random brightness variations significantly larger than the measurement uncertainty. Second, the typical peak brightness of the robustly determined variables is shifted to the right (brighter), compared with the full sample even when only considering protostars (red dots). This is an expected occurrence due to the increased signal-to-noise ratio for bright sources and its carryover effect on the completeness of our analyses. We return to this discussion in Section 4.

In this section, we consider the quantifiable properties of the robustly recovered variables in our sample in order to determine an approximate classification scheme. Using a false alarm threshold criterion of 0.1% and a stochastic threshold of ± 5σ, we recover only 18 secular variables among the 295 bright submillimeter sources in our sample.

All the robust variables within our survey are known protostars. None of the known disk sources show any variability in the submillimeter using our LSP, linear, and stochasticity analyses. This absence is at least in part due to the low peak brightness of the disk sources in submillimeter (see Figure 7), and the brightness-dependent completeness discussed in Section 4. We further discuss the variability depending on evolutionary stage in Section 5.2. Limiting our sample only to protostars, we find that ∼22% (18 out of 83) are observed to be secular variables and none show evidence of a single stochastic outlier event.

Tables 3 and 4 present the quantitative results for individual secular variables. The observed amplitude (ΔF) is defined as F_max–F_min where F_max and F_min are the brightest and faintest peak brightness among the measured fluxes of a source. The LSP-derived periods of the detected variables vary from 1 year to 40 years. We categorize the variables with their best-fit periods compared to the 4 yr observing window of our survey: Periodic (<4 yr: Group P), Curved (4–15 yr: Group C), and Linear (>15 yr: Group L).

Table 3. Physical Properties of Robust Secular Variables

Region	ID	Figure a	Known Name	$\bar{{\rm{F}}}$	ΔF/σfid	ΔF/$\bar{{\rm{F}}}$	Lbol	Tbol
				(Jy bm−1)		(%)	(L⊙)	(K)
IC 348	J034356.5 + 320050	B.1	IC348 MMS 1	1.39	4.57	10.2	1.6	23 b
Serpens main	J182951.2 + 011638	B.2	EC 53	1.17	16.2	37.7	5.9	161 b
NGC 2068	J054647.4 + 000028	B.3	HOPS 389	0.99	5.09	12.4	6.0	43 c
NGC 1333	J032910.4 + 311331	B.4	IRAS 4A	9.16	7.84	15.7	8.3	31 b
Serpens main	J182952.0 + 011550	B.5	Serpens SMM 10	0.84	6.21	16.1	8.3	70 b
Ophiuchus	J162626.8-242431	B.6	VLA 1623-243	3.89	5.42	11.0	0.5	45 d
NGC 2068	J054613.2-000602	B.7	V1647 Ori	0.25	7.93	47.3	27	322 c
IC 348	J034357.0 + 320305	B.8	HH 211	1.21	6.59	15.2	1.4	23 b
NGC 2068	J054603.6-001447	B.9	HOPS 315	0.52	3.06	10.2	6.2	180 c
NGC 2068	J054631.0-000232	B.10	HOPS 373	1.26	10.9	25.0	5.3	37 c
OMC2/3	J053529.8-045944	B.11	HOPS 383	0.54	5.56	18.1	7.8	46 c
NGC 1333	J032903.8 + 311449	B.12	West 40	0.48	7.86	27.8	0.7	18 b
Serpens south	J182937.8-015103	B.13	IRAS 18270-0153	0.63	5.79	17.2	6.9	57 b
Serpens main	J182949.8 + 011520	B.14	Serpens SMM 1	7.08	9.35	18.8	70	13 b
Serpens main	J182948.2 + 011644	B.15	SH 2-68 N	2.10	3.20	6.76	14	31 b
Serpens south	J183004.0-020306	B.16	CARMA 7	4.85	5.31	10.7	15	35 e
NGC 2068	J054607.2-001332	B.17	HOPS 358	1.29	16.2	37	25	42 c
NGC 1333	J032903.4 + 311558	B.18	NGC1333 VLA 3	3.04	11.3	23	38	220 b

Notes.

^a The reference number of the light curve in Figure Set B (online-only material). ^b For these sources, the L_bol and T_bol values are taken from Mowat (2018). See also Section 5.1 for more details. ^c For these sources, the L_bol and T_bol values are taken from Furlan et al. (2016). ^d For this source, the L_bol and T_bol values are taken from Murillo et al. (2018). ^e For this source, the L_bol and T_bol values are taken from Maury et al. (2011).

Download table as:  ASCIITypeset image

Table 4. Light-curve Properties of Robust Secular Variables

Known Name	Period a	A/σfid	A/$\bar{{\rm{F}}}$	Slope/$\bar{{\rm{F}}}$	Grp b	log FAPMod	log FAPLin
	(yr)		(%)	(% yr−1)
MMS 1	1.05	1.78	2.70	−1.3	P	−3.36	−3.27
EC 53	1.61	8.77	10.8		P	−10.7	−0.29
HOPS 389	4.53	2.69	3.17	−1.5	C	−4.90	−5.39
IRAS 4A	4.54	0.46	4.1	−2.0	C	−3.79	−3.04
SMM 10	5.30	4.50	4.61		C	−4.48	−0.40
VLA 1623-243	5.72	0.668	2.57	−1.6	C	−3.19	−3.72
V1647 Ori	5.83	30.0	18.1	−10.4	C	−14.0	−12.2
HH 211	5.84	2.86	3.81	2.15	C	−7.56	−7.30
HOPS 315	8.16	3.79	3.26		C	−3.14	−2.32
HOPS 373	8.16	8.94	12.2		C	−7.27	−0.75
HOPS 383	8.16	9.56	8.11	−2.80	C	−8.08	−5.39
West 40	8.18	11.0	8.55	4.84	C	−8.37	−9.15
IRAS 18270-0153	12.5	12.1	11.1	−2.0	C	−4.56	−3.68
SMM 1	37.1			2.07	L	−5.91	−4.69
SH 2-68 N	37.1			−0.96	L	−4.22	−3.85
CARMA 7	37.5			1.29	L	−4.99	−5.36
HOPS 358	40.8			−7.06	L	−12.5	−13.7
VLA 3	40.9			3.34	L	−7.10	−7.79

Notes.

^a The best-fit period from the LSP method. ^b P: Periodic Group. C: Curved Group. L: Linear Group.

Download table as:  ASCIITypeset image

We find two protostars in the Periodic Group: the known periodic infrared and submillimeter variable source EC 53 in Serpens Main (Hodapp et al. 2012; Yoo et al. 2017; Lee et al. 2020b) with ∼1.5 yr period, and the protostellar source MMS1 in Perseus IC 348 with ∼1 yr period. MMS1 is an intriguing source because it is also detected through our linear analysis (see the light curve presented in Figure B1.1). The majority, 61%, of the recovered secular variables belong to the Curved Group, for which the observing window is comparable to the derived source periodicity. As mentioned earlier, for these sources, the quantified period is more likely a measure of variability timescale rather than a determination of a sinusoidal nature. These sources require dedicated, long-term monitoring over many years to further unravel their episodic nature. Finally, the Linear Group contains five protostars, 28% of the secular variables, that show predominantly a linear trend in their light curves but with the possibility of slight curvature. We discuss each of these individual sources in the Appendix.

Alternatively to studying periodicity timescales, we can also consider the distribution of the linear fits to the light curves to determine if sources are in a brightening or dimming phase. The distribution of slopes, in units of fractional change per year, is obtained from the least-squares linear fitting to the light curves of all protostars and is shown in Figure 6 (light red) overlaid with the 14 robust linear fits (red hatched). For the protostars with robust fits, nine are declining in brightness while only five are rising. Furthermore, the typical strength of the decline or brightening is similar, although there are two significant outliers with large negative slopes. The lack of symmetry between rising and dimming submillimeter sources is intriguing but suffers from small number statistics at present.

Zoom In Zoom Out Reset image size

Across the four years of observations analyzed here, the distribution of the slopes over all investigated protostars (light red) is strongly peaked near the origin, and has a mean value of 0.0016, essentially zero, and σ_obs(4 yr ) = 0.013, equivalent to a rising or falling slope of 1.3% per year. Individually, these slopes are not considered robust precisely because the measured slope is smaller than, or similar to, the uncertainty in the slope determination. The ensemble, however, contains useful information on the allowable range of these individually uncertain slopes. Previously, on the basis of the first 18 months of observations, Johnstone et al. (2018) fit a Gaussian distribution to a similarly derived set of light-curve slopes and measured σ_obs(1.5 yr ) = 0.023, of which at least half of the spread was calculated to be due to the uncertainty in fitting the slopes, an uncertainty that decreases with longer observing timelines. Recognizing that measurement uncertainty contributed to their observed spread in slopes, Johnstone et al. (2018) determined that the underlying intrinsic distribution of light-curve slopes is σ_int ≲ 0.01 and predicted that, after three years of observations, the observed width of the distribution should shrink to be at most σ_obs(3 yr ) ∼ 0.01, with the upper limit occurring if σ_int = 0.01. Here, when considering all the protostars in our sample, we obtain that upper limit. We therefore predict that submillimeter monitoring surveys that can reach a relative brightness calibration significantly below 1% will find that a majority of the protostellar sources are robust secular variables.

Johnstone et al. (2018) recovered five robust secular variables within the JCMT Transient Survey sample after the first 18 months of the survey. All of those sources continue to be identified as robust secular variables by the present study. Furthermore, the lengthening of the observing window to 4 yr has increased the number of recovered secular variables by greater than a factor of three.

Finally, as mentioned above, the robust secular variables are biased toward the brighter protostars. In Figure 7, we present histograms of the peak brightness of our full source sample as well as subsets for protostars, disks, and secular variables. We discuss the issue of completeness within our secular variability analysis in the next section.

Zoom In Zoom Out Reset image size

Although our survey reveals only a small number of variables, which limits the statistical significance of results from subsamples, we consider possible regional variations in the numbers of secular variables uncovered. Table 6 shows the ratio (${P}_{\sec }$) between the numbers of secular variables and protostars in each region, to examine any regional dependency. For the full survey, we find 22% of protostars varying, while across the eight individual star-forming regions, we find a wider range: 0%–50%. To account for the completeness (Section 4), we repeat our analysis using only those sources with mean peak brightness greater than 0.5 Jy beam⁻¹ and find that 37% of protostars vary across the entire survey.

We compare the regional variability ratios by calculating the probability (p) that each region is drawn from an underlying sample with the mean value. Each region is compared against the entire sample, excluding the region of interest. When taking the samples to include all protostars, there is no clear evidence of regional dependence. Serpens Main is somewhat overabundant in variables and has p = 0.04 (P < 0.05 is typically used as a dissimilarity threshold), but has also been most frequently observed and should have a slightly higher completeness threshold than the other regions. Considering only protostars brighter than 0.5 Jy beam⁻¹, where the completeness is more uniform, the OMC 2/3 region becomes somewhat exceptional in harboring fractionally few variables, 1 out of 12 protostars, with p = 0.01. There is no obvious correlation in these numbers with distance to each region. The full set of results is presented in Table 6.

Protostars are expected to be more variable at earlier evolutionary stages. For this analysis, we use the bolometric temperature, corrected for extinction, as a proxy for evolutionary stage. For the three Orion fields, these values were directly obtained from the literature (Furlan et al. 2016). For Serpens, Perseus, and Ophiuchus, bolometric temperature and luminosity values were updated by Mowat (2018) based on the SED-fitting methodology in Dunham et al. (2015) using SCUBA-2 fluxes at 450 and 850 μm from the JCMT Gould Belt survey (Pattle et al. 2015; Chen et al. 2016) to improve the submillimeter coverage. ³⁷

In Figure 10, we plot both the extinction-corrected bolometric luminosities (upper panel) and the peak submillimeter brightnesses (lower panel) of our sample against the extinction-corrected bolometric temperature. We also include the larger Orion HOPS sample (Furlan et al. 2016) in the luminosity panel.

Zoom In Zoom Out Reset image size

Considering just the luminosity versus temperature (upper panel), the protostellar variables are more luminous and cooler than the underlying HOPS sample, but not significantly more luminous than the typical protostar in the JCMT Transient Survey sample. However, in the bottom panel it is clear that peak submillimeter brightness is a strong function of evolutionary class, as expected. Care must be taken before inferring changes in the variability fraction with evolutionary stage due to changing completeness limits as a function of source brightness.

Following the results of Section 4, we consider only those sources for which the mean brightness is greater than 0.5 Jy beam⁻¹, for which the completeness to secular variability is reasonably flat. We thus find 0 (out of 2) Class II secular variables, 3 (out of 10; 30%) Class I secular variables, and 13 (out of 34; 38%) Class 0 secular variables. Although the variability ratio is higher among Class 0 than Class I, the difference between the ratios of Class 0 and Class I is insignificant (p = 0.61).

We conducted a literature search for known eruptive variables within the eight JCMT Transient Survey fields. Six such stars were identified across four regions, of which three are coincident with submillimeter peaks. All three of the coincident submillimeter protostars, EC 53 in Serpens Main (also known as V371 Serpentis; Hodapp et al. 2012), V1647 Ori in NGC 2068 (Aspin & Reipurth 2009), and HOPS 383 in OMC 2/3 (Safron et al. 2015) are observed to be varying (see Table 4). Based on bolometric temperature (see Section 5.2), the first two are classified as a Class I while HOPS 383 is a Class 0. The three known eruptive stars for which there is no coincident submillimeter peak are OO Serpentis ³⁸ and V370 Serpentis, both in Serpens Main, as well as IRAS 18270-0153W ³⁹ in Serpens South.

The two optical/IR eruptive (FUOr/EXOr) stars that are found to vary in our submillimeter sample, V1647 Ori and EC 53, are both observed to have episodic accretion outbursts. EC 53 has a short-term, ∼1.5 yr, periodicity, observed at near-IR and submillimeter wavelengths (see, Lee et al. 2020b). V1647 Ori has undergone multiple eruptions, with times between bursts in the 2 to 5 year range (Aspin et al. 2006; Acosta-Pulido et al. 2007; Ninan et al. 2013). Most recently, in 2018, V1647 Ori was observed in a quiescent phase (Giannini et al. 2018). The decrease during 2018 is observed in our survey (see Figure B1.3). The submillimeter light curve for this source is fit with a ∼6 yr period, placing it in the Curved Group. The light curve itself suggests that it may be reaching a minimum (Figure B1.3).

These two known optical/IR eruptive Class I sources are among the most extreme submillimeter variables in our sample in terms of fractional flux change across the observing window. Only HOPS 383, our mid-infrared (MIR) identified eruptive Class 0 variable, along with the three Class 0 sources HOPS 358, HOPS 373, and West 40, reveal a similarly large brightness range (see also Section 5.6). Furthermore, similar to HOPS 383, the three strongly varying Class 0 sources have extrapolated timescales (periods) greater than 8 yr. Interestingly, HOPS 383 has a somewhat lower variability amplitude across the last 4 years, compared with the other three Class 0 variables. This may be due to the fact that its submillimeter brightness decay has recently slowed, after the well-known outburst event between 2004 and 2010 (Safron et al. 2015).

Additional details on each of these protostellar variables can be found in the Appendix.

For embedded protostars, their luminosity is dominated by accretion luminosity. Therefore, the low-luminosity sources, which are classified as YSOs with ≤1 L_⊙ (Dunham et al. 2008), should have low mass accretion rates. The most extreme low-luminosity sources, called Very Low-Luminosity Objects (VeLLOs), were discovered by the Spitzer survey (Young et al. 2004). By definition, VeLLOs are embedded protostars with the internal luminosity ≤0.1 L_⊙, and are thus considered as YSOs in the most quiescent phase of the episodic accretion process. We investigate the variability of low-luminosity sources and VeLLOs revealed in submillimeter.

First, we cross-matched those low-luminosity sources from Dunham et al. (2008) with the source list in this study. We note that the SMM 1 was classified as a low-luminosity source by Dunham et al. (2008) using the 2MASS and Spitzer data set (1.25–70 μm). ⁴⁰

There are 28 low-luminosity sources located within our coverage, with seven LLSs brighter than >0.14 Jy beam⁻¹ at 850 μm: two in IC 348, four in NGC 1333, and one in Ophiuchus. The two low-luminosity sources in IC 348, MMS 1 and HH 211, are secular variables (Figures B1.1 and B1.8). Thus, two out of seven low-luminosity sources (∼30%) show secular variability in the submillimeter. Comparing this variability fraction against the entire sample of monitored protostars (see Table 6), we obtain a p-value of 0.69. For the sources F₈0% > 0.5 Jy beam⁻¹, the number of variables from the sample becomes two out of four (50%), for which the p-value is 0.59. Thus, we do not detect any clear differentiation in the secular variability fraction of low-luminosity sources, with any conclusions limited by the small sample size.

As an alternate sample, we compare our source list to the VeLLO list from Kim et al. (2016). Kim et al. (2016) added four complementary criteria (a ratio of the 1.65 μm flux to the 70 μm flux less than 2.8; a 4.5 μm magnitude brighter than 15.3; not being registered as galaxies in known databases; and a color index [8]–[24] > 2.2) from Dunham et al. (2008) for identifying VeLLOs. Only four VeLLOs from this sample, all in the Serpens South region, are detected in the submillimeter, although 14 VeLLOs (one in IC 348, two in NGC 1333, two in Ophiuchus, and nine in Serpens South) are within our coverage. None of the matched VeLLOs shows secular variability in our analysis. Despite the null detection of variability among these four observed VeLLOs, there is no evidence that the small VeLLO sample is different than the larger survey sample. That is, comparing the fractional variability within the VeLLO sample against the larger survey sample yields a p-value ≈0.14, well above the 0.05 typically used to justify differences in the samples.

According to our results, submillimeter continuum emission traces only variability in protostellar luminosity, which appears as a temperature change in the thick envelope. Optical and near-IR brightness variations are sensitive to the luminosity, but may also be caused by extinction changes within the small inner region close to the central object, (i.e., optical depth to the stellar photosphere and inner disk). As a corollary, the variability of YSOs at later evolutionary stages, without much remaining envelope, is primarily revealed at shorter wavelength. Furthermore, different physical and geometrical scales associated with underlying variability processes also determine the timescales of variability.

In this sense, the MIR covers a wider range of physical/geometrical scales, and thus, a greater variety of variability mechanisms than are expected to be traced in the submillimeter. An intensive investigation for YSO variability in the MIR has been undertaken for ∼5400 YSOs in 20 nearby low-mass star-forming regions, using NEOWISE W2 band (4.6 μm) light curves covering ∼7 years (Park et al. submitted). Here, we compare the JCMT Transient Survey results to that MIR variability investigation.

In total, 1204 NEOWISE sources are located within our eight Transient survey regions. Of those, 49 have submillimeter counterparts: 39 are classified as protostars, 9 are disk sources, and 1 submillimeter peak is coincident with a Class III or evolved source. For the protostars in common, 23 out of 39 (59%) are variables in the MIR over the 6.5 years monitored by NEOWISE, very similar to the 55% variability likelihood over all protostars detected by that survey (W. Park et al. submitted).

From our survey, 8 submillimeter protostellar variables have counterparts observed in the W2 band, yielding a 20% (8 out of 39) variability rate almost identical to the likelihood over the full submillimeter survey. Five of these sources are identified as variable at both wavelengths. West 40, V1647 Ori, and Serpens Main-SMM 1 are classified as Curved in both submillimeter and MIR. EC 53 and Serpens Main-SMM 10 are varying in W2 but classified as MIR irregulars, which means no notable secular trend in the W2 light curves. HOPS 389, NGC 1333-VLA 3, and SH 2-68-N are not observed to vary in the MIR.

Similar to the submillimeter analysis in this work, Park et al. (submitted) classified the YSO MIR variability largely into two categories, i.e., secular and stochastic. As here, they further divided secular variability into three types: Linear, Curved, and Periodic; however, the boundary between the groups is different because of the different cadences and coverage in the two surveys. Given the large number of MIR stochastic variables that they found, they further divided the stochastic variability into three types: Burst, Drop, and Irregular. Burst and Drop are identified by sudden brightening and dimming only in a few epochs (i.e., with short timescales) over the 7 yr light curve, while Irregular is identified by the random distribution of brightness with a high standard deviation, in which no underlying timescale of variability is measured.

The distinct difference between protostellar variability observed in the MIR versus the submillimeter is that most of the MIR variability is irregular while all variables in the submillimeter show secular trends. This difference can be interpreted in terms of both the different cadences of observations and the different origins for the variability at these two wavelengths. The cadence of NEOWISE is about 6 months, but that of the JCMT Transient Survey is typically one month or shorter. As a result, periodic MIR variability with a short period is more likely to be classified as irregular variability for the NEOWISE light curves. EC 53 provides a good example; it is classified as an irregular variable in the MIR, but as a periodic variable in the submillimeter. The period of EC 53 is about 1.5 yr. However, the periodogram analysis in the MIR (W. Park et al. submitted) did not find it as a periodic source because the light curve is not simply sinusoidal (see Figure 1 and Lee et al. 2020b), and the MIR data cadence is too sparse to pull out the periodicity.

There is an additional explanation for the significantly larger numbers of observed irregular variables in the MIR compared with the submillimeter. MIR brightness is intrinsically sensitive to the variability of the protostar, while the submillimeter radiation is sensitive only to the variability of the envelope temperature, where the submillimeter emission arises. As a result, stochastic changes in the accretion luminosity or extinction events taking place close to the central source can be more easily detected in the MIR. Sudden changes in accretion luminosity, however, are smoothed by the large envelope, given its long thermal relaxation time of a month or longer (Johnstone et al. 2013), and the fractional change in submillimeter flux, which is proportional to the temperature change, is lowered by a factor of ∼5.5 compared to that in MIR (Contreras Peña et al. 2020).

Finally, there are no epochs in which the observed submillimeter emission becomes much lower, i.e., drops, for a single epoch. In the MIR sample, a non-negligible number of disk sources show such short-timescale behavior. The lack of dips in the submillimeter is most likely due to the fact that the observed submillimeter radiation emitted by the envelope traces the time-averaged luminosity of the source, where the averaging is over the months of thermal equilibrium time of the radiating envelope (see Johnstone et al. 2013). Thus, while it is possible that a nonthermal brightening event such as a flare (Mairs et al. 2019) might add to this emission to produce a burst, there are no short-timescale subtractions of emission. Furthermore, the submillimeter radiation from the envelope is optically thin, and thus responds only linearly to changes in optical depth, unlike the optical through MIR where small changes in the line-of-sight column density can provide very large optical depth variations and strong dimming on a variety of timescales.

A key result from four years of monitoring protostars in the submillimeter is that at least 20% are seen to undergo secular variations in their brightness. The fraction rises to 40% when only the brighter protostars are included in the sample. Nevertheless, the measured fractional variability amplitudes and slopes for individual sources (see Table 4, columns 6 and 7) appear to be modest, and thus one might be tempted to dismiss the observed variability as physically unimportant with respect to the mass assembly process for protostars. However, as noted by Johnstone et al. (2013), the submillimeter emission reacts to the temperature of the protostellar envelope, effectively probing the Rayleigh–Jeans tail of the Planck curve, and thus converting the observed variability to the underlying protostellar brightness variations requires a significant, exponential boost (see also MacFarlane et al. 2019a, 2019b). Through comparison of observed variables at both MIR and submillimeter wavelengths, Contreras Peña et al. (2020) empirically uncovered a submillimeter-to-protostellar luminosity exponential factor of ∼4, and confirmed that it matched well with radiative transfer expectations (Baek et al. 2020). Thus, assuming that the protostellar luminosity is a reasonably linear proxy for protostellar accretion, it is appropriate to consider the importance of the variable mass accretion uncovered by these submillimeter observations.

We begin by assuming that the 12 secular variables best fit with periodic light curves undergo bursts with duration timescales half the measured period and burst amplitudes equal to twice the amplitude, A_burst = 2 × A (see Table 4), effectively requiring that quiescence is determined at the lowest point on the observed light curve. Furthermore, for those six variables best fit by linear slope measurements, we assume that the burst takes place during half of the observing window, ∼2 yr, and that the burst amplitude is four times the (yearly) slope, A_burst = 4 × Slope (see Table 4). This latter assumption is likely a significant underestimate of the importance of the burst, as we are constraining the amplitude by the observing time window.

With these estimates for the submillimeter burst strength, we can now directly compute the excess amount of mass accreted during each burst after subtraction of the quiescent mass accretion, M_burst, against the total mass that is accreted at the fixed quiescent rate over the full time period, M_quiescent:

$$\begin{eqnarray}&&\displaystyle \frac{{M}_{\mathrm{burst}}}{{M}_{\mathrm{quiescent}}}=0.5({\left[1+{A}_{\mathrm{burst}}\right]}^{4}-1).\end{eqnarray} \tag{ 8 }$$

Thus, a null value refers to no variability, A_burst = 0, whereas a very powerful burst would add many multiples of the quiescent accreted mass. For the quantities tabulated in Table 4, we find that the smallest amplitude variables add ∼10% of the quiescent accreted mass (integrated over the full time range) through their years-long burst enhancement. The six most prominent variables, however, each add greater than 40% during the burst phase.

Using this simple model, the most powerful burst accretor in our sample is V1647 Ori, adding an additional 125% of the quiescently accreted mass during each ∼3 yr burst. Furthermore, during each burst V1647 Ori attains a peak mass accretion rate 3.5 times the quiescent rate,

$$\begin{eqnarray}&&\displaystyle \frac{{\dot{M}}_{\mathrm{peak}}}{{\dot{M}}_{\mathrm{quiescent}}}={\left[1+{A}_{\mathrm{burst}}\right]}^{4}.\end{eqnarray} \tag{ 9 }$$

Similarly, EC 53 (V371 Ser) adds an additional 60% of the quiescently accreted mass during each of its ∼ 0.75 yr bursts, attaining a peak burst mass accretion rate twice the quiescent rate. As noted in Section 5.3, both of these known Class I repeaters belong to the FUOr/EXOr type of eruptive stars and have repetition timescales in the literature similar to what is determined in our analysis.

In addition to EC 53 and V1647 Ori, four Class 0 protostars in our sample have measured high amplitude variability and associated large burst accretion rates: HOPS 356 (85%), HOPS 373 (70%), West 40 (45%), and HOPS 383 (40%) (see Table 4). All of these sources are associated with long timescales, periods of at least twice the 4 yr monitoring window, and thus the measured amplitudes and burst accretion rates are anticipated to be lower limits.

The simple analysis presented in this section suggests that the mass assembly of protostars can be significantly influenced by bursts on timescales of years. Roughly 7%, 6 out of the 83 monitored protostars, show order unity mass accretion variability on this timescale. Assuming that all protostars sample the range of variability observed by our sample, we predict that, for a given source, years-long bursts occur on average about every ∼50 yr. Thus, the importance of these accretion events for the overall mass assembly of the typical protostar is not large, accounting for only a few percent of the accreted mass.

The results presented here are in contrast to the extreme bursting events searched for in the MIR by Scholz et al. (2013) and Fischer et al. (2019), where less than 1% of the protostars were seen to vary on timescales of half a decade, implying an underlying random distribution of extreme bursts for a given source with ∼1000 yr between each event. There remains a vast unexplored time domain separating these two extremes for episodic mass accretion that warrants detailed investigation. This will require continued searches for rare brightening events as well as dedicated monitoring of protostellar samples in order to determine if the variability uncovered in this paper smoothly connects with the rarer but more extreme events or if the variability separates cleanly into well-defined episodic timescales. Both possibilities would yield significant constraints for theoretical models of mass accretion, both steady-state and unstable, through the disk.

Finally, we note that continued monitoring in the submillimeter is essential for the youngest, most embedded, protostars. The four strongly varying Class 0 sources identified here, with variability amplitudes similar to the known eruptive systems EC 53 and V1647 Ori, are hard to interpret from MIR observations alone, due to the extreme extinction of the central source and disk (see, e.g., the multiwavelength analysis of HOPS 373 by Yoon et al. in prep.). ALMA observations are also revealing complicated geometries during these very early phases of protostar evolution, such as multiple protostars and an arc-like structure within an extremely dense and compact core (Tokuda et al. 2014) and multiple misaligned outflows from a single young protostar (Okoda et al. 2021). The explanation for these features is assumed to be complex gas accretion within a turbulent core. Such complex structure will inevitably produce protostellar accretion variability. It will, therefore, be important to also study the relation between submillimeter variability and the small-scale structure near the protostar using high spatial resolution observations such as by ALMA.

The submillimeter light curve for secular source Perseus IC 348, MMS 1 is presented in Figure B1.1. The JCMT source is associated with a millimeter peak with a separation of 066 from the Class 0 protostar IC 348 MMS. Continuum observations by by Chen et al. (2013) with the SMA established IC 348 MMS as a multiple system, with two continuum sources, MMS1 and MMS2, embedded in a common envelope and separated by 98. IC 348 MMS1 drives molecular outflow traced in H₂ (Eislöffel et al. 2003) with distinct red and blueshifted lobes detected in CO (Lee et al. 2016). Additionally, Segura-Cox et al. (2016) identified MMS1 as a disk candidate in the VLA Nascent Disk and Multiplicity (VANDAM) Survey. The second continuum source, MMS2, is potentially a protobinary confirmed by presence of an outflow, HH797, with signatures of a rotating jet (McCaughrean et al. 1994; Pech et al. 2012).

The submillimeter light curve for secular source Serpens Main EC 53 (V371 Serp) is presented in Figure B1.2. Based on its bolometric temperature and spectral index, EC 53 is a Class I source (Evans et al. 2009; Dunham et al. 2015). However, the envelope contains ∼6 M_⊙ (Baek et al. 2020), while the protostar and disk masses are only 0.3 M_⊙ and 0.07 M_⊙, respectively (Lee et al. 2020a), suggestive of a much younger evolutionary state. Periodic variability has been observed previously in the near-IR (Hodapp 1999; Hodapp et al. 2012), MIR (Contreras Peña et al. 2020), and the submillimeter (Mairs et al. 2017b; Yoo et al. 2017; Johnstone et al. 2018). Deep CO and H₂O absorption lines, indicative of viscous disk heating and characteristic of FUOrs, are also observed (Connelley & Reipurth 2018). Lee et al. (2020b) combine these observations across wavelengths to analyze the physical process responsible for each burst.

The submillimeter light curve for secular source Orion NGC 2068, HOPS 389 is presented in Figure B1.3. Launhardt et al. (1996) observed this source in the 1.3 mm continuum using the IRAM 30 m telescope, finding it to have a compact condensation with more extended envelope. The presence of two embedded infrared sources led to its identification as a self-luminous protostellar core (Class 0/I). Mitchell et al. (2001) mapped the source using SCUBA and weakly detected an associated outflow in ¹²CO 3–2. Walker-Smith et al. (2013) also observed patches of red- and blueshifted ¹³CO 3–2 emission in the vicinity of this source. Matthews & Wilson (2002), using the SCUBA polarimeter, found the polarization toward the source to be consistent with that of its surroundings.

Spezzi et al. (2015) list one source for these coordinates, classified as Class I based on the SED slope between 2.2 and 24 μm. Furlan et al. (2016) find two closely associated sources, HOPS 323 and 389, classifying HOPS 323 as Class I and HOPS 389 as Class 0. Here, we identify the submillimeter source with HOPS 389, although HOPS 323 is somewhat more luminous. Both sources are listed as protostars in the Spitzer survey of YSOs in Orion (Megeath et al. 2012). These two objects, together with HOPS 322, are described as a group of protostars.

The VLA/ALMA Nascent Disk and Multiplicity (VANDAM) Survey of Orion Protostars (Tobin et al. 2020) detected four localized 870 μm sources, presumed to be disks: HOPS-323-A at 05:46:47.697 + 00:00:25.27 (44.475 ± 0.754 mJy) and HOPS-323-B at 05:46:47.667 + 00:00:24.81 (81.351 ± 1.484 mJy), both Class 1; and HOPS-389-A at 05:46:47.019 + 00:00:27.07 (19.722 ± 0.668 mJy) and HOPS-389-B at 05:46:46.604 + 00:00:29.52 (4.785 ± 0.325 mJy), both Class 0.

The submillimeter light curve for secular source Perseus NGC 1333, IRAS 4A is presented in Figure B1.4. The bright infrared (Dunham et al. 2015) through submillimeter (Kirk et al. 2006; Enoch et al. 2009) source is classified as Class 0. IRAS 4A is resolved as a double source by Tobin et al. (2016a) and appears to harbor a massive disk ∼1 M_⊙ (Tychoniec et al. 2018) and a large-scale bipolar molecular outflow (Choi 2001). The source was monitored by Spitzer at 3.6 and 4.5 μm for ∼35 days as part of YSOVAR (Rebull et al. 2015), but was not found to be variable in the MIR.

The submillimeter light curve for secular source Serpens Main, SMM 10 is presented in Figure B1.5. The source has been identified as Class I (Kryukova et al. 2012), and is also known as Ser-emb-12 (Enoch et al. 2009). Near-infrared K-band observations have identified a stellar counterpart to the submillimeter source, SMM 10 IR, which is found to vary by about 2 mag in K band, and is associated with a fan-shaped nebulosity (Hodapp 1999). Correlated secular variability on year-to-year timescales is found for both MIR NEOWISE and Transient Survey submillimeter observations of the source (Contreras Peña et al. 2020). ALMA 12 m array band 6 observations (Project ID: 2013.1.00618.S; ALMA source name: 211) of the dust continuum at 094 resolution reveal the source to be a binary, with a separation of 425, hosting a complex outflow structure in the ¹²CO line.

The submillimeter light curve for secular source Ophiuchus Core, VLA 1623-243 is presented in Figure B1.6. VLA 1623-243 comprises three protostars: VLA 1623 A, an embedded Class 0 source 56 from the peak of 850 μm emission; VLA 1623 B, a Class 0 source 65 from the peak; and VLA 1623 W, a Class I YSO that is much farther away, with projected distance of 160 (Ward-Thompson et al. 2011; Murillo & Lai 2013). VLA 1623 A and B are just 12 apart, and thus cannot be resolved by the JCMT. We therefore associate the binary system with the observed submillimeter light curve, although VLA 1623 A is brighter (Murillo et al. 2018) and thus expected to contribute most to the observed flux. VLA 1623 A has an envelope mass ∼0.8 M_⊙ (Murillo et al. 2018), a disk with an outer radius of 180 au traced by C¹⁸O emission (Murillo et al. 2013), and a large-scale bipolar outflow (Santangelo et al. 2015). VLA 1623 B is a weaker source with about a quarter of the envelope mass of that surrounding VLA 1623 A (Murillo et al. 2018). Santangelo et al. (2015) also found a smaller strongly collimated outflow associated this component. Interestingly, VLA 1623 B was determined to be variable at centimeter wavelengths by Coutens et al. (2019).

The submillimeter light curve for secular source Orion NGC 2068, V1647 Ori is presented in Figure B1. V1647 Ori is the only viable target that is close to the observed submillimeter peak. It is classified as FUOr-like (Connelley & Reipurth 2018) with multiple bursts: 2003, 2008, and 2011 (Aspin et al. 2006; Acosta-Pulido et al. 2007; Ninan et al. 2013; Giannini et al. 2018). The source was reported to dim in July 2012 by Ninan et al. (2012), four years before the start of the JCMT Transient Survey. Its classification, based on the spectral energy distribution, changes over time, but recent work (Furlan et al. 2016) classifies it as a Class I source in transition to Class II with a flat spectrum. Andrews et al. (2004) reported an increase of a factor of 25 at 12 microns and an order-of-magnitude increase in bolometric luminosity up to 34 L_⊙. In the quiescent phase, Ábrahám et al. (2004) approximate V1647 Ori's luminosity around 5–6 times solar.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (18 panels)

Tarball of figure set images

Using VLTI, Ábrahám et al. (2006) found no companion around V1647 Ori. There is a possible flared disk around V1647 Ori with a disk mass of 0.05 M_⊙ as inferred from MIR. Mosoni et al. (2013) modeled the MIDI observations toward this source and found that the disk+envelope pre- and post- burst structures are similar. With the recent ALMA observations, Cieza et al. (2018) constrain the target to a 40 au disk whose mass is 80 Jupiter masses and with an inclination of 57°. They did not include the envelope model in this case. The earliest jets observations were obtained by Eisloffel & Mundt (1997). A molecular CO outflow was later detected by Andrews et al. (2004). From infrared imaging, the outflow direction is north to south with respect to the source.

The long-term evolution of this source up to 2006 is reported by Acosta-Pulido et al. (2007). The I-band magnitude dropped by 5 magnitude from 2004 February to 2006 January, with a rapid drop at the end of 2005. A periodicity of 56 days was determined from the long-term decay that can be explained by dust structures as indicated by Eiroa et al. (2002). A high-cadence photometry reported by Bastien et al. (2011) uncovered flickers that may be due to interactions between the magnetosphere and the disk. Aspin (2011) reported a new flare in 2011 with a total luminosity of 16 L_⊙ and an accretion rate of 4 × 10⁻⁶ M_⊙ yr⁻¹, which is similar to the earlier flare in 2004 (Fedele et al. 2007).

The submillimeter light curve for secular source Perseus IC 348, HH 211 is presented in Figure B1.8. The JCMT peak coincides with the dense core HH 211. This continuum source comprises two respective embedded sources, HH 211 SMM1 and HH 211 SMM2, with an angular separation of 031 resolved by high-resolution observations with the SMA (Lee et al. 2009; Chen et al. 2013).

HH 211 SMM1 is an embedded young Class 0 protostar object, potentially the youngest in Perseus (Hirota et al. 2008). It is the exciting source (Avila et al. 2001) for a highly collimated jet observed in H₂ (McCaughrean et al. 1994) and SiO (Hirano et al. 2006). Furthermore, Gueth & Guilloteau (1999) detected molecular outflow activity in CO, with low-velocity CO emission tracing bipolar cavities and high-velocity emission aligned with the collimated jet-like structure closer to the central object. Lee (2020) notes that the knots of emission in the jet have estimated timescales between decades and centuries. The central source appears as a flattened disk-like structure, with a size ∼80 au, inclined perpendicular to the jet and outflow axis (Gueth & Guilloteau 1999; Lee et al. 2009).

Wiseman (2001) uncovered a rotating envelope in ammonia emission, normal to the outflow and very close to the central protostar. Additionally, Tanner & Arce (2011) and Tobin et al. (2011) detected notable line broadening indicative of outflow–envelope interaction using ammonia observations. High-resolution mapping with ALMA, by Lee et al. (2018), confirmed the presence of a disk and found small outflow velocities indicating a rotating disk atmosphere. Properties of the adjacent continuum source, SMM2, which has a lower mass than SMM1, are mostly uncertain.

The submillimeter light curve for secular source Orion NGC 2068, HOPS 315 is presented in Figure B1.9. HOPS 315 is classified as a Class I protostar via SED fitting by Furlan et al. (2016). Könyves et al. (1919) identified a dense core 36 away from our JCMT peak using PACS and SPIRE images from the Herschel Gould Belt survey. The authors derived a critical Bonnor–Ebert mass ratio of 0.2 for the core, indicating that it is self-gravitating. Recently, Tobin et al. (2020) used ALMA 870 μm observations to derive protostellar disk properties, including the dust disk radius ∼46 au and dust disk mass ∼100 Earth masses. Kounkel et al. (2016) found no companions within 100–1000au of HOPS 315 using near-IR observations. ALMA observations of the source and outflow are also presented by Dutta et al. (2020).

The submillimeter light curve for secular source Orion NGC 2068, HOPS 373 is presented in Figure B1.10. Gibb & Little (2000) identified this source as a Class 0 protostar, and mapped a bipolar outflow with high-velocity dense gas. The geometric center of this outflow was close to the position of a previously detected H20 maser (Haschick et al. 1983). Edwards & Snell (1984) had previously detected high-velocity molecular gas in the vicinity of the maser. Phillips et al. (2001) observed the source with SCUBA and found its 450 μm emission to be elongated perpendicular to the outflow direction. In CARMA 2.9 mm continuum imaging (Tobin et al. 2015), the source appears to be a binary. The outflow was further studied with CARMA (Tobin et al. 2016b), and a compact second component with the opposite orientation was seen, perhaps due to a second outflow from the binary source. From the (Herschel PACS) O I (63.18 μm) line intensity, they estimate a mass-loss rate of 1.1 × 10⁻⁷ M_⊙ yr⁻¹.

The VLA/ALMA Nascent Disk and Multiplicity (VANDAM) Survey of Orion Protostars (Tobin et al. 2020) detected two 870 μm sources, presumed to be disks: HOPS-373-A at 05:46:31.099-00:02:33.02 (103.827 ± 1.626 mJy) and HOPS-373-B at 05:46:30.905-00:02:35.19 (92.206 ± 0.947 mJy), classifying both as Class 0.

The JCMT Gould Belt Survey (Kirk et al. 2016) did not find an associated protostar in the Spitzer (Megeath et al. 2012) or Herschel (Stutz et al. 2013) survey catalogs of YSOs in Orion. However, the source is included in the second paper on the Spitzer Orion survey (Megeath et al. 2016) as a newly Spitzer-identified YSO (there classified as disk) and later listed by Stutz et al. (2013) as a PACS Bright Red source (PBRS).

The source was observed using the SCUBA polarimeter by Matthews & Wilson (2002), and while polarization was seen along the filament where it is located, no polarization was detected toward the source itself. The authors suggest that this end of the filament may lie in the foreground of the nebula, as do Walker-Smith et al. (2013).

Kang et al. (2015) find a high deuteration ratio (HDCO/H2CO) for this source, implying that it is probably in the very earliest stage of star formation. Its classification as a PBRS also argues for it being a very young protostar, and this is corroborated by the models of Furlan et al. (2016). However, Tobin et al. (2015) found rapidly declining visibility amplitudes and suggest it may be slightly more evolved.

The submillimeter light curve for secular source Orion Orion A OMC2/3, HOPS 383 is presented in Figure B1.11. Safron et al. (2015) reported an apparent MIR outburst of HOPS 383 between 2004 and 2006. By 2008, its brightness at 24 μm became 35 times brighter than the observation in 2004, and subsequent monitoring suggested no significant decrease in luminosity from 2009 to 2012. In 2017, Fischer & Hillenbrand (2017) suggested a decline in the NIR luminosity. The authors used SED modeling to predict the NIR flux (J, H, and K_s band), and found no detection at HOPS 383 in the NIR imaging. Based on the post-outburst SED, HOPS 383 is classified as a Class 0 protostar (Safron et al. 2015).

Galván-Madrid et al. (2015) found no corresponding radio outburst between 1998 and 2014 after searching the VLA archives for radio counterparts of HOPS 383. In late 2017, Grosso et al. (2020) detected a hard X-ray source with Chandra, which had not been observed in an earlier, 2000 January, epoch. The newly found X-ray source is spatially coincident with the radio source (JVLA-NW) imaged by Galván-Madrid et al. (2015). Meanwhile, Grosso et al. (2020) monitored the same region at near-IR wavelengths and identified an H₂ knot ∼15'' away. The H₂ knot was cross-matched with a previously observed source, SMZ 1-2B, obtaining a proper motion of 18 in the southeastern direction over 20 years. According to this proper motion, Grosso et al. (2020) estimated the dynamical timescale of the outflow shocked H₂ knot at 180 ± 100 years.

The submillimeter light curve for secular source Perseus NGC 1333, West 40 is presented in Figure B1.12. Also referred to as Per-emb-15 (Enoch et al. 2009), West 40 has been observed from the infrared (Dunham et al. 2015) through the radio, where it is unresolved by the VLA (Tobin et al. 2016a). The source is classified by Dunham et al. (2015) as Class 0.

The submillimeter light curve for secular source Serpens South, IRAS 18270-0153 is presented in Figure B1.13. The location of the submillimeter peak has not been observed by either the SMA or ALMA. The object is identified as a Class 0 source by Dunham et al. (2015). As noted in the footnote to Section 5.3, Johnstone et al. (2018) previously misidentified this source as the FU Ori candidate IRAS 18270-0153W (Connelley & Greene 2010), but closer inspection shows that it is beyond 15'' from that infrared-bright source.

The submillimeter light curve for secular source Serpens Main, SMM 1 is presented in Figure B1.14. SMM 1 (Casali et al. 1993) is also known as Serpens FIRS 1 (Harvey et al. 1984) and Ser-emb-6 (Enoch et al. 2009). This source is the brightest Class 0 in the Serpens Main region and one of the most extensively studied. Observations with the 12 m ALMA array reveal two extremely bright resolved sources, SMM1-a and the relatively fainter SMM1-b, surrounded by complex extended structure associated with outflow cavities (Hull et al. 2016; Francis et al. 2019). Hull et al. (2016) find high-velocity ∼80 km s⁻¹ CO jets emanating from SMM1-a and -b, and interpret a C-shaped structure in the dust continuum around SMM1-a as walls of a cavity carved by precession of the jet. The same cavity is also seen in free–free emission in VLA observations, which Hull et al. (2016) suggest to be caused by ionization of gas in shocks at the cavity walls. Polarization measurements with ALMA suggest that the jets are also playing a role in shaping the local magnetic field (Hull et al. 2017). Lower-velocity (∼10–20 km s⁻¹) wide-angle outflows are also seen in the CO emission around the high-velocity jets (Hull et al. 2014, 2017). Studies of the chemistry of the outflows have identified SiO in the high-velocity jet and wide-angle outflows, while HCN and H₂CO are only seen in the slower wings, consistent with a lower C/O ratio in the jet (Tychoniec et al. 2019). Mid-infrared Spitzer observations also find jets in H₂ and various atomic emission lines (e.g., [Fe ii]); however, interpreting which source is driving each outflow is complicated by the complexity of the outflows and the lower Spitzer resolution (Dionatos et al. 2014).

The submillimeter light curve for secular source Serpens Main, SH 2-68 N is presented in Figure B1.15. SH 2-68 N, also known as Ser-emb 8 (Enoch et al. 2009) is an embedded protostar within an extended structure that also includes Ser-emb 8N. This Class 0 + 1 source lies about 2'' north of the bright submillimeter peak (Francis et al. 2019). Hull et al. (2017) reported on ALMA observations of polarized dust emission. They found a weak and randomly oriented magnetic field, on scales of 100–1000 au, which did not exhibit the expected hourglass shape, likely due to the enhanced role of turbulence at these large scales. Sh 2-68 N has a slow bipolar molecular outflow observed in SiO 5–4 (Hull et al. 2014) powered by an extremely high-velocity jet. The outflow has a wide opening angle in ¹²CO 2–1, surrounded by SO emission tracing the cavity walls (see also Aso et al. 2019).

The submillimeter light curve for secular source Serpens South, CARMA 7 is presented in Figure B1.16. The submillimeter peak lies in a crowded region, right in the middle of the Serpens South Cluster. The brightest source in the JCMT field seen with ALMA is serp45 (Plunkett et al. 2018), and we identify the JCMT source with it. This source is also known as CARMA 7 (3 mm; Plunkett et al. 2015a) and VLA 12 (Kern et al. 2016). This source is a Class 0 protostar with clear evidence for an episodic jet (Plunkett et al. 2015b). The source was not identified as a protostar by Dunham et al. (2015); however, Herschel observations reveal that the source is both luminous and buried within a cold envelope (identified as SerpS-mms18 by Maury et al. 2011).

The submillimeter light curve for secular source Orion NGC 2068, HOPS 358 is presented in Figure B1.17. The variable submillimeter source is associated with the protostar HOPS 358. Stutz et al. (2013) classified this source as a PACS Bright Red sources (PBRS) with a 70 to 24 μm flux ratio greater than 1.65. The authors suggest that PBRs sources are extreme Class 0 objects with a higher envelope density compared with typical Class 0 sources, and equivalently, higher rates of mass infall. Furlan et al. (2016) also classified this source as a Class 0 protostar via SED fitting. Nagy et al. (2020) investigated the outflow properties, finding a dynamical timescale of 10⁴ yr and a mass-loss rate of 5 × 10⁻⁶ M_⊙ yr⁻¹. The authors also showed an infall asymmetry. ALMA observations of the source and outflow are also presented by Dutta et al. (2020). Tobin et al. (2020) derived protostellar dust disk properties, finding a roughly 135 au dust disk.

The submillimeter light curve for secular source Perseus NGC 1333, VLA 3 is presented in Figure B1.18. The submillimeter source was previously identified by Enoch et al. (2009) as Per-emb-44. It is a Class I protostar according to Dunham et al. (2015). Although identified as a single source by Tobin et al. (2016a), it separated into two peaks at cm wavelengths wavelengths (Tychoniec et al. 2018). The source was monitored by Spitzer at 3.6 and 4.5 μm for ∼35 days as part of YSOVAR (Rebull et al. 2015), but was not found to be variable in the MIR on that timescale.