The JCMT Transient Survey: Identifying Submillimeter Continuum Variability over Several Year Timescales Using Archival JCMT Gould Belt Survey Observations

The accretion history of low-mass stars has been under investigation for many years (see, e.g., Kenyon et al. 1990; Bell & Lin 1994). Despite recent advances (Vorobyov & Basu 2005; Rice et al. 2010; McKee & Offner 2011; Dunham & Vorobyov 2012; Bae et al. 2014; Vorobyov & Basu 2015), few observational constraints exist regarding the dominant mass transfer processes occurring during the earliest phases of star formation. The presence of a disk around a protostar complicates how matter is moved from the nascent envelope to the central source. Through a combination of rotation and magnetic fields, anisotropies develop in the otherwise symmetric infalling material predicted in the seminal model of Shu (1977), causing this material to first accrete onto the disk before being transported to the protostar (for a review on pre-main sequence accretion, see Hartmann et al. 2016). The rate of this mass transport likely varies on both long and short timescales, due to gravitational (Vorobyov & Basu 2005, 2006), magnetorotational (Armitage et al. 2001; Zhu et al. 2009), and spiral wave (Bae et al. 2016) instabilities. Material that builds up in the disk can compress the magnetosphere of a forming star, leading gas to rapidly flow onto the star through directed funnels (Romanova et al. 2003; see also Cody et al. 2017). In addition, the formation of giant planets (Nayakshin & Lodato 2012) and stellar encounters in binary systems (Forgan & Rice 2010; Hodapp et al. 2012) can also affect the rate at which a protostar gains mass.

Episodic accretion (short bursts of mass accretion separated by long, quiescent periods) caused by these instabilities in the disk has been gaining considerable attention in the literature (see Audard et al. 2014, and references therein). These events may be related to known outbursts observed in FU Orionis (Herbig 1977; Hartmann & Kenyon 1985; Reipurth 1990) and EX Lupi (Herbig 1989) objects. In addition to these large outbursts, Classical T Tauri Stars (CTTS) have been known to exhibit lower levels of variability (e.g., Johns & Basri 1995; Alencar et al. 2001; Bouvier et al. 2007; Donati et al. 2013; Blinova et al. 2016). Accretion variability events with analogous intensities and timescales may occur during the youngest stages of a forming star's evolution. A non-steady accretion rate is also one solution to "the luminosity problem": the order of magnitude discrepancy between the fainter, observed median protostellar luminosity and the predicted protostellar luminosity assuming constant, steady state accretion (Kenyon et al. 1990). Other solutions to this problem, such as longer timescales over which low-mass stars form, have also been posed (e.g., McKee & Offner 2011). A recent study of 3000 Young Stellar Objects (YSOs) present within 18 different molecular clouds conducted by Dunham et al. (2015) using data obtained by the Spitzer Space Telescope, however, has provided evidence that Class 0+I protostellar lifetimes are relatively short (0.46–0.72 Myr; see also Evans et al. 2009; Heiderman & Evans 2015; Ribas et al. 2015; Fischer et al. 2017).

Observing the changes in accretion rate around a forming star is fundamental to understanding the earliest stages of star formation. This is possible by measuring significant changes in the protostar's peak brightness over time, since the protostar's luminosity is generated primarily by the accreting material. Johnstone et al. (2013) modelled the spectral energy distribution (SED) of a deeply embedded protostar undergoing an outburst due to an increased mass accretion rate and concluded that while a stronger signal will be detected at the peak of the SED (mid- and far-infrared wavelengths), the outburst also would be observable in the submillimeter regime. The submillimeter luminosity arises from cascade reprocessing of the stellar radiation by the surrounding dust. The negligible heat capacity for interstellar dust causes the time lag between the burst and the observation at submillimeter wavelengths to be dominated by the light crossing time of the protostellar envelope, which Johnstone et al. (2013) calculated to be several weeks to months long.

The benefit of observing variability at submillimeter wavelengths is that ground-based instruments such as the James Clerk Maxwell Telescope (JCMT) are available to monitor large areas of the sky at a regular cadence. The JCMT Transient Survey (Herczeg et al. 2017) is obtaining simultaneous continuum data at 450 and 850 μm on eight nearby (<500 pc) star-forming molecular clouds at approximately 28 day intervals, when the targets are observable using the Submillimeter Common-User Bolometre Array 2 (SCUBA-2). These eight 30' in diameter regions were originally selected from the JCMT Gould Belt Survey (GBS; Ward-Thompson et al. 2007) and shifted on sky to optimize the number of Class 0+I and Class II YSOs present in each field. The GBS data that overlap with the Transient Survey data were collected between 2012 and 2014, while the Transient Survey data collection began in 2015 December. Therefore these two surveys can be compared to identify and investigate variable signals over 2–4 year timescales. Here we consider Transient Survey data obtained prior to 2017 March 1. In the present work, we only use 850 μm data, which has higher signal to noise ratio (S/N) than the 450 μm data (see Section 2).

The GBS data were not obtained with the goal of detecting protostellar variability, so the observations are not regularly spaced in time, and they often occur over only one or two nights for a given field. Therefore in this study we measure the flux-calibrated peak brightnesses of extracted sources in the co-added GBS images and compare them with the same sources in the co-added Transient Survey images. In this way, we robustly characterize the properties of all identified sources and become more sensitive to long-term changes. Furthering our investigation, we also correlate the identified 850 μm emission sources with the positions of known YSOs (Megeath et al. 2012; Stutz et al. 2013; Dunham et al. 2015). This is the first study that systematically analyses two consistent submillimeter data sets separated in time in order to determine constraints on the magnitude and the frequency of deeply embedded, variable protostars.

This paper is organized as follows: In Section 2, we summarize the Transient Survey, GBS, and YSO catalog observations we use throughout this work. In Section 3, we give an overview of our data reduction, image alignment, and relative flux calibration procedures. In Section 4, we present the results of the comparison between the GBS and Transient Survey. This includes the identification of variable candidates, a description of the quality of those candidates, and the calculated fractional peak brightness changes per year. In Section 5, we discuss the results, highlight previously known variable sources in our sample, and construct a light curve for the variable source EC 53 (Hodapp et al. 2012; Yoo et al. 2017). Finally, in Section 6, we summarize our main results and provide concluding remarks.

In this analysis, we combine JCMT Transient Survey (Herczeg et al. 2017) observations with archival data of the same regions observed by the JCMT Gould Belt Survey (hereafter GBS; Ward-Thompson et al. 2007). All observations in both surveys are performed using the Submillimeter Common-User Bolometer Array 2 (SCUBA-2) instrument (Holland et al. 2013) at the JCMT. SCUBA-2 provides continuum coverage at both 450 μm and 850 μm simultaneously, with half-power bandwidths of 32 μm and 85 μm (Holland et al. 2013) and effective beam sizes of 98 and 146, respectively (Dempsey et al. 2013). While we expect a stronger variability signal toward the peak of the SED, typically the mid to far-IR (Johnstone et al. 2013), the practical data reduction and analysis of SCUBA-2 450 μm images, is complicated by a strong atmospheric water vapor dependency, in addition to a less stable beam profile due to dish deformation and focus errors (Dempsey et al. 2013). These effects dramatically influence the S/N and require further investigation before careful corrections are possible. In this paper we focus only on the 850 μm data. All of the observations were taken in the PONG1800 mapping mode (Kackley et al. 2010), yielding circular maps ("PONGs") ∼05 in diameter.

In total, there are eight Transient Survey fields (Herczeg et al. 2017): three in the Orion A and B Molecular Clouds (OMC 2–3, NGC 2068, and NGC 2024), two in the Perseus Molecular Cloud (IC 348 and NGC 1333), two in the Serpens Molecular Cloud (Serpens Main and Serpens South), and one in the Ophiuchus Molecular Cloud (Oph Core). Eleven Gould Belt Survey fields include bright, compact sources within areas of significant overlap of submillimeter emission with the Transient Survey data. Central locations, observation dates, 850 μm noise, and program identification numbers for all observations are summarized in Table 1 (see also, Appendix A).

Table 1.  A Summary of the Observed Fields and Their Co-added Noise at 850 μm

Tile Namea	Central R.A. (J2000)	Central Decl. (J2000)	Start Dateb	End Dateb	${\rm{\Delta }}t$ c	850 μm Noised (Jy beam−1)	# Obs.e	SCUBA-2 Program Identification
IC348	03:44:18.3	32:05:16.0	2015 Dec 22	2017 Feb 09	⋯	0.0043	9	M16AL001
IC348-E	03:44:24.4	32:02:08.6	2012 Aug 16	2013 Feb 05	3.5	0.0055	4	MJLSG38
NGC 1333	03:28:54.5	31:17:09.0	2015 Dec 22	2017 Feb 06	⋯	0.0039	10	M16AL001
NGC1333-N	03:29:06.3	31:22:26.7	2012 Jul 02	2012 Jul 03	4.0	0.0052	4	MJLSG38
OMC2–3	5:35:33.2	−5:00:32	2015 Dec 26	2017 Feb 06	⋯	0.0042	9	M16AL001
OMC1 tile4	5:35:49.4	−4:46:23	2012 Aug 17	2013 Aug 26	3.2	0.0045	4	MJLSG31
NGC 2068	5:46:13.0	−0:06:05	2015 Dec 26	2017 Feb 06	⋯	0.0039	10	M16AL001
OrionBN_450_S	5:46:13.0	−0:06:05	2014 Nov 16	2014 Nov 22	1.6	0.0046	6	MJLSG41
NGC 2024	5:41:41.0	−1:53:51	2015 Dec 26	2017 Feb 06	⋯	0.0043	11	M16AL001
OrionBS_450_E	5:42:48.0	−1:54:36	2014 Oct 27	2014 Nov 09	1.6	0.0045	6	MJLSG41
OrionBS_450_W	5:40:32.2	−1:48:00	2013 Feb 12	2013 Mar 03	3.3	0.0045	4	MJLSG41
Oph Core	16:27:03.2	−24:32:46.5	2016 Jan 15	2017 Feb 06	⋯	0.0050	8	M16AL001
L1688-1	16:27:02.9	−24:41:44.5	2012 May 06	2012 May 08	4.0	0.0053	4	MJLSG32
L1688-2	16:27:15.1	−24:10:09.7	2012 May 18	2012 May 20	4.0	0.0051	4	MJLSG32
Serpens Main	18:29:48.7	01:15:39.5	2016 Feb 02	2017 Feb 22	⋯	0.0045	9	M16AL001
SerpensMain1	18:29:59.8	01:14:46.9	2012 May 18	2012 May 19	4.1	0.0058	4	MJLSG33
Serpens South	18:30:02.2	−02:02:23.0	2016 Feb 02	2017 Feb 22	⋯	0.0047	9	M16AL001
SerpensS-NE	18:31:35.4	−01:53:50.3	2012 Apr 21	2012 May 03	4.2	0.0049	4	MJLSG33
SerpensS-NW	18:29:30.8	−01:47:07.3	2012 May 03	2012 May 05	4.2	0.0048	5	MJLSG33

Notes. JCMT Transient Survey fields are in bold, and associated JCMT Gould Belt Survey fields are listed below each Transient Survey field.

^aIn nearly all cases, the tile name corresponds to the target identifier in the JCMT archive. The exception is NGC2068, which is labeled as NGC 2071 in the archive. ^bThe start and end dates refer to the date of the first observation taken by each survey and the last observation taken before 2017 March 1 (yyyymmdd). ^cThe time between the average GBS date and the average Transient Survey date (years). ^dThese measurements of the 850 μm noise levels are based on a point source detection in each field's co-added image after smoothing with a 6'' FWHM Gaussian kernel. The effective beam size after smoothing is 158. ^eThe number of observations included in the co-add.

Download table as:  ASCIITypeset image

2.1. The JCMT Transient Survey

2.2. The JCMT Gould Belt Survey

2.3. Spitzer Space Telescope and Herschel Space Observatory YSO Catalogs

We largely follow the three-step data reduction and calibration methodology adopted by the JCMT Transient Survey, which is described in detail by Mairs et al. (2017). We develop further methods, however, described in Section 3.2, to perform a relative flux calibration between the two independent data sets (the GBS and the Transient Survey). We first construct robust images from the raw SCUBA-2 data for all observations. Then we perform a spatial alignment for each image using a reference field. Finally, we perform a relative flux calibration to bring all observations into agreement with the mean, co-added image (for the GBS and Transient Surveys separately).

3.1. Data Reduction

3.2. Post-reduction Alignment and Flux Calibration

Nominally, the JCMT has a pointing error of 2''–6'' and a flux calibration uncertainty of ∼5%–10% (Dempsey et al. 2013; Mairs et al. 2017). Using the methods developed by Mairs et al. (2017), we achieve an alignment of $\lt 1^{\prime\prime}$ and a flux calibration uncertainty (in a typical measurement) of 2%–3%. Since each image has a similar flux calibration uncertainty, co-adding the data further reduces the uncertainty in a given measurement by the square root of the number of observations included in the co-add. Therefore a source's peak brightness may have an uncertainty as low as 1% or less (see Tables 1 and 2). Briefly, the image alignment and flux calibration procedures rely on comparing the properties of bright, compact sources detected across multiple observations of the same area of the sky. To identify these sources; subtract larger-scale background flux; and extract properties such as the central position, size, and peak brightness, we employ the algorithm Gaussclumps (Stutzki & Guesten 1990), which models each compact object with a Gaussian profile. Specifically, we use the starlink software (Currie et al. 2014) implementation of Gaussclumps found within the cupid (Berry et al. 2007) package. We expect the bright, point-like sources in an image to resemble Gaussian structures based on the shape of the JCMT beam. In order to mitigate spurious noise features, before running the source extraction algorithm, we first smooth the maps with a 6'' (FWHM) Gaussian kernel (see Appendix B of Mairs et al. 2017).

Table 2.  Summary of the Variable Candidate Source Properties

GBS Field	ID	Other Namea	R.A. (J2000)	Decl. (J2000)	${f}_{\mathrm{trans}}$ b	$\tfrac{{\sigma }_{{f}_{\mathrm{trans}}}}{{f}_{\mathrm{trans}}}$ c	${f}_{\mathrm{GBS}}$ b	$\tfrac{{\sigma }_{{f}_{\mathrm{GBS}}}}{{f}_{\mathrm{GBS}}}$ c	$| \delta |$	$\tfrac{\dot{f}}{{f}_{\mathrm{trans}}}$ (% yr−1)	${\sigma }_{\dot{f}/f}$	Category
NGC1333-N	PER-1	IRAS4Ad	3:29:10.42	31:13:30.63	8.83	1.0	8.10	0.4	7.66	2.09	0.28	S
NGC1333-N	PER-10	Bolo 40e	3:28:59.86	31:21:33.09	0.60	1.0	0.67	0.7	7.99	−3.03	0.36	S
NGC1333-N	PER-34	[LAL96] 213f	3:29:07.66	31:21:54.05	0.28	2.8	0.38	1.2	8.31	−8.64	0.87	S
SerpensMain1	SER-1	[KOB2004] 258bg	18:29:49.80	1:15:19.33	6.38	1.5	5.76	0.4	6.85	2.37	0.39	S
SerpensS-NW	SER-10	IRAS 18270-0153h	18:29:37.99	−1:51:04.66	0.66	0.8	0.78	0.6	11.81	−4.05	0.35	S
L1688-2	OPH-14	⋯	16:26:24.95	−24:24:23.92	0.42	1.9	0.32	2.7	8.04	6.10	0.78	E
SerpensS-NW	SER-21	SerpS-MM15i	18:30:02.69	−2:01:09.33	0.42	0.9	0.37	1.3	7.32	2.82	0.38	E
OMC1 tile4	ORA-36	HOPS 383j	5:35:29.67	−4:59:37.25	0.53	0.9	0.58	1.6	4.17	−2.66	0.64	Pos.

Notes.

^aReference name from the literature. YSO name where possible. ^b ${f}_{{\rm{x}}}$ is the mean source peak brightness measured across the Transient Survey (x = trans) or GBS (x = GBS) in Jy beam⁻¹. ^cThe standard deviation in f divided by the square root of the number of observations, normalized by f. In units of %. ^dJennings et al. (1987). ^eEnoch et al. (2006). ^fLada et al. (1996). ^gKaas et al. (2004). ^hConnelley et al. (2007). ⁱMaury et al. (2011). ^jSafron et al. (2015).

Download table as:  ASCIITypeset image

As described by Mairs et al. (2017), the image alignment is performed by selecting the robust, Gaussian sources that have a peak brightness of at least 200 mJy beam⁻¹ (S/N ∼ 10 for an individual Transient Survey observation) and a maximum effective radius of 10'' (the effective radius is $r\,=\sqrt{{\mathrm{FWHM}}_{1}\times {\mathrm{FWHM}}_{2}}/2$, where the FWHM_N terms are the full widths at half maximum of the major and minor axes of the fitted two-dimensional Gaussian), comparing their central locations in each observation to a reference image, and then correcting for that offset. For both the Gould Belt and Transient Surveys, the chosen reference field for each region is the first Transient Survey observation of that region. The absolute position of each source is of little importance for our goals, as the most critical measurement for understanding the variability of a given point source is the relative peak brightness and the alignment uncertainty is small enough that we are able to confidently associate known protostars with the emission peaks.

Mairs et al. (2017) describe the procedure to identify and use robust calibrator sources in each field to self-consistently perform a relative flux calibration. The chosen calibrator sources are referred to as Family Members. A given Family Member in the Transient Survey is selected by measuring the peak brightness of all sources that are brighter than 500 mJy beam⁻¹ with effective radii $\lt 10^{\prime\prime}$ normalized to their average peak brightness and comparing that value with all of the other sources in a given image with these properties. The largest set of sources that display a low amount of scatter from observation to observation with respect to one another (defined by a threshold of 6% in the standard deviation) is selected to be a Family. The fact that these sources agree well with one another over time suggests that none of them are intrinsically varying to the level we can detect and that they are tracking the flux uncertainty of the telescope. In this manner, we determine if each observation, as a whole, is slightly brighter or fainter than the mean image allowing each epoch (image) to be corrected by a constant multiplicative factor.

We perform the relative flux calibration individually for the GBS data and the Transient data. As described in Section 2, the GBS data were often taken over a short time frame (one or two days), whereas the Transient data were taken over the course of more than 1 year. Therefore sources that do not appear to vary across the GBS observations could show signs of variability in the Transient data. In addition, a slowly varying source may look constant in both data sets individually, but given the long separation in time between the two surveys (see Table 1), the peak brightness could be dramatically different when cross-compared (see Section 4).

In order to bring the GBS data and the Transient Survey data into relative flux calibration with one another, we first identify sources common to both surveys. The sources we select have peak brightnesses larger than $200\,\mathrm{mJy}\ {\mathrm{beam}}^{-1}$ and effective radii less than $10^{\prime\prime}$. For each source, we measure the average peak brightness (${f}_{\mathrm{GBS}}$ and ${f}_{\mathrm{trans}}$) and derive the associated uncertainties (${\sigma }_{{f}_{\mathrm{GBS}}}$ and ${\sigma }_{{f}_{\mathrm{trans}}}$) by calculating the measured standard deviation across all observations within each survey individually and dividing the result by the square root of the number of observations that were included in the calculation (see Table 2).

We calculate a calibration factor using the measured source brightnesses by dividing the ${f}_{\mathrm{GBS}}$ values by their corresponding ${f}_{\mathrm{trans}}$ values. For the ith source, we label this ratio R_i and we propagate the uncertainties, ${\sigma }_{R,i}$, in order to calculate the weighted mean, $\bar{R}$,

$$\begin{eqnarray}&&\bar{R}=\displaystyle \frac{{\displaystyle \sum }_{i=1}^{n}{R}_{i}{w}_{i}}{{{\rm{\Sigma }}}_{i=1}^{n}{w}_{i}},\end{eqnarray} \tag{ 1 }$$

where n is the number of sources and the weights are the inverse square of the propagated uncertainties,

$$\begin{eqnarray}&&{w}_{i}={\sigma }_{R,i}^{-2}.\end{eqnarray} \tag{ 2 }$$

The uncertainty in the weighted mean, ${\sigma }_{\bar{R}}$, is given by

$$\begin{eqnarray}&&{\sigma }_{\bar{R}}=\sqrt{\displaystyle \frac{1}{{{\rm{\Sigma }}}_{i=1}^{n}{w}_{i}}}.\end{eqnarray} \tag{ 3 }$$

Next, we calculate the difference, δ, between each peak brightness ratio and the weighted mean by subtracting the latter from the former and dividing by the uncertainty associated with that source added in quadrature to the uncertainty in the weighted mean,

$$\begin{eqnarray}&&{\sigma }_{\mathrm{tot},i}=\sqrt{{\sigma }_{R,i}^{2}+{\sigma }_{\bar{R}}^{2}},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{\delta }_{i}=\displaystyle \frac{{R}_{i}-\bar{R}}{{\sigma }_{\mathrm{tot},i}}.\end{eqnarray} \tag{ 5 }$$

We consider all sources with $| \delta | \gt 3$ to be deviant outliers that might skew the weighted mean. Therefore we remove these sources from the calculation and recompute the weighted mean. This refined weighted mean (that excludes outliers) represents the initial approximation of the relative flux calibration factor between overlapping average GBS and Transient Survey images.

To explore how the peak brightness uncertainties in each data set (${\sigma }_{{f}_{\mathrm{GBS}}}$ and ${\sigma }_{{f}_{\mathrm{trans}}}$) affect the refined weighted mean, we use a Monte Carlo analysis. We fix ${\sigma }_{{f}_{\mathrm{GBS}}}$ and ${\sigma }_{{f}_{\mathrm{trans}}}$ for each source and randomly draw new peak brightness measurements from normal distributions with mean values equal to ${f}_{\mathrm{GBS}}$ and ${f}_{\mathrm{trans}}$ and standard deviations equal to ${\sigma }_{{f}_{\mathrm{GBS}}}$ and ${\sigma }_{{f}_{\mathrm{trans}}}$. In this way, we calculate a new peak brightness ratio for each source, R_i, that is within the derived measurement uncertainties. Then, we calculate a new weighted mean (Equation (1)) based on these values, discard outliers in the same manner as before, and compute a refined weighted mean. We repeat this process 10,000 times.

Of these 10,000 refined weighted mean values, we adopt the average as the relative flux calibration factor, or relative FCF. This is the number by which we divide the GBS image to bring it into relative calibration with the Transient Survey image. The standard deviation in refined weighted mean values is the uncertainty in the relative FCF, ${\sigma }_{\mathrm{FCF}}$.

We plot the results in Figures 1, 2, and 3. The colors represent a source's association with known YSOs. In these figures, to be associated with a protostar, flat spectrum source, protostellar candidate, or disk, the peak position of the source must be within 10'' of the YSO location (to match the radius of the largest compact source we consider). Points labeled with a "c" are Family Members in both the GBS and Transient Surveys (i.e., they appear to have stable peak brightnesses in each data set independently but not necessarily when compared over several year timescales). The FCF values are indicated by the dashed lines, and ${\sigma }_{\mathrm{FCF}}$ is indicated by the dotted lines. As the noise is relatively constant across a map and large-scale modes have been mostly removed in the data reduction procedure, we expect the relative FCF value to remain constant across the full field of any single image.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The GBS images were all originally calibrated at a level that is slightly brighter than their respective Transient Survey images obtained prior to 2017 March 1.This is the result of several factors that affect the nominal flux calibration performed at the JCMT, before we apply the relative flux calibration presented previously. For instance, the nominal flux calibration values are based on data obtained between 2011 and 2012 (Dempsey et al. 2013), the water vapor monitor was replaced in 2015, and SCUBA-2 had a filter upgrade in late 2016. The combination of these factors has led to differences in the original calibration of the data obtained in the GBS and the Transient Survey eras. These effects are corrected for by our relative flux calibration.

In Section 3.2, we described how we bring the mean GBS data and the mean Transient Survey data into relative flux calibration with one another by using the bright, compact sources in each overlapping field. Once the relative FCF has been computed, we repeat our calculation to search for significant outliers by replacing the individually calculated weighted mean of the source peak ratios, $\bar{R}$, with the mean FCF and ${\sigma }_{\bar{R}}$ with ${\sigma }_{\mathrm{FCF}}$ in Equation (5). Here we define the significance of outliers based on the distribution of these newly calculated δ values.

Figures 4, 5, and 6 show the δ value of every source. In Figure 7, we show the distribution of δ values for all sources fit with a Gaussian curve. The fitted Gaussian has a standard deviation of ${\sigma }_{\mathrm{GaussFit}}=1.76$. The measured δ values significantly deviate from the fitted Gaussian profile beyond a threshold of $\delta =\pm 6$. Therefore we define sources that have $| \delta | \gt 6$ as noteworthy outliers. The probability of detecting a source with $| \delta | \gt 6$ from the Gaussian distribution of values is 0.06%. Over all 11 GBS fields that overlap with the Transient Survey data, however, we find that 7 out of 175 independent sources²⁰ brighter than 200 mJy beam⁻¹ with radii <10'' exceed this threshold. Thus 96% of the identified sources show no sign of variation between the average GBS data and the average Transient data to our sensitivity, while the remaining 4% are considered variable candidates and form the basis of all further investigation. We summarize the properties of the variable candidates in Table 2.

In order to discern the quality of each of the variable candidates, we construct difference maps, subtracting the GBS co-add from the Transient co-add for each region, and perform a visual analysis on each map to note sources with extended features in confused regions as well as larger-scale differences between the maps which can complicate our measurements (see Appendix A).

While investigating the difference maps, we searched for any indication of compact sources that were present only in the GBS data and not in the Transient Survey data (or vice versa), in case an object had varied such that it fell below the detection threshold (200 mJy beam⁻¹) in one data set but not in the other. No significant objects of this type were identified.

Through our visual inspection of each source in their respective difference maps, we define 2 categories of variable candidate:

• 1.  
  Strong Candidates (S): The source exceeds the significance threshold ($| \delta | \gt 6$), has an average Transient Survey peak brightness measurement of ${f}_{\mathrm{trans}}\gt 200\,\mathrm{mJy}\ {\mathrm{beam}}^{-1}$, and has a radius less than 10''. These sources have an obvious indication of significant compact structure at the location of the source in the difference map.
• 2.  
  Extended Candidates (E): The source exceeds the significance threshold ($| \delta | \gt 6$), has an average Transient Survey peak brightness measurement of ${f}_{\mathrm{trans}}\gt 200\,\mathrm{mJy}\ {\mathrm{beam}}^{-1}$, and a radius less than 10''. There is little indication of compact structure present at the location of the source in the difference map; surrounding extended structures complicate the measurements.

We find 5 Strong sources and 2 Extended sources. Extended variable candidates display very little evidence of real compact emission changing brightness between the GBS and Transient Survey data sets. It is likely that we are underestimating the uncertainty for these objects due to the uncertainty in the background subtraction computed by the Gaussclumps algorithm in confused areas. Both of these remain interesting sources, though they are less robust candidates for variability than those classified as Strong. The Extended source OPH-14 traces part of the larger-scale structure around the prototypical Class 0 source VLA 1623 (Andre et al. 1993), though we see no significant evidence that the deeply embedded protostar itself is undergoing a significant brightness change at 850 μm. In addition to the Strong and Extended candidates, one Possible candidate (HOPS 383; Safron et al. 2015) is also identified in Tables 2 and 3, and is discussed further in Section 5.1.

In general, we expect the presence of a disk to influence the accretion rate of material onto the central protostar. Of the five Strong variable candidates, four are associated with a known protostar (and likely, therefore, a young circumstellar disk), and none are associated with a known evolved disk object (see N_proto in Table 3; see also Table 4). To be associated with a protostar or a disk, the peak position of the source must be within 10'' of the YSO location. Recall that the presence of a dusty envelope increases the likelihood of detecting variability at submillimeter wavelengths due to the reprocessing of the light emitted by a burst event. Generally the envelope is faint or non-existent for known, more evolved disk candidates seen at infrared wavelengths.

In total, out of the 175 independent sources identified across all 8 Transient Survey fields and their 11 associated GBS fields, there are 54 protostars and 13 disk objects within 10'' of source peaks. Therefore, to the sensitivity achieved across these surveys, approximately 7% of known protostars associated with an 850 μm emission source are shown to be potentially varying. Only 3% of the sources not known to be associated with a YSO show signs of variability. One of the Strong variable candidates, IRAS4A (PER-1), has two protostars associated with its peak. Another, SER-1, displays an elongated structure that may indicate multiple embedded sources or some outflow activity. Excluding the Extended sources, which are less robust detections, only one source without a known protostar or disk shows signs of variability. We do not expect starless cores to be variable. Therefore the Strong variable candidate without a known embedded YSO, Bolo 40 (Enoch et al. 2006), is a good target for follow-up studies to identify signatures of very faint, deeply embedded protostars that have yet to be detected.

The GBS data and the Transient Survey data were obtained 2–4 years apart, depending on the field (see Table 1), which allows us to characterize the apparent, significant submillimeter brightness changes in terms of a fractional change per year (averaged over both data sets), assuming a linear change throughout the GBS/Transient Survey time lag,

$$\begin{eqnarray}&&\displaystyle \frac{\dot{f}}{{f}_{\mathrm{trans}}}=\displaystyle \frac{{f}_{\mathrm{trans}}-{f}_{\mathrm{GBS}}}{{f}_{\mathrm{trans}}}\times \left(\displaystyle \frac{1}{{\bar{t}}_{\mathrm{trans}}-{\bar{t}}_{\mathrm{GBS}}}\right),\end{eqnarray} \tag{ 6 }$$

where ${\bar{t}}_{\mathrm{GBS}}$ is the mean GBS observation date for that field and ${\bar{t}}_{\mathrm{trans}}$ is the mean Transient Survey observation date for that field. (The parenthetical term in the denominator is the same as ${\rm{\Delta }}t$ in Table 1.) We choose to normalize the result to the mean Transient Survey observation, to express the brightness change as a percentage of the (near) current source peak brightness. The associated error, ${\sigma }_{\dot{f}/f}$, is the combination of the Transient Survey average peak brightness error (see Section 3.2 and Table 2), the GBS peak brightness error, and the error in the weighted mean.

We plot $\dot{f}/{f}_{\mathrm{trans}}$ as a function of the Transient Survey average peak brightness in Figure 8, colored according to category. A typical, Strong variable candidate has an average brightness change per year of $| 4.0| \% \,{\mathrm{yr}}^{-1}$, with a standard deviation of $2.7 \% \,{\mathrm{yr}}^{-1}$. The Extended variable candidates are grouped in a similar area of parameter space. This grouping is primarily the result of poor background flux subtraction performed by the Gaussclumps algorithm in these more confused regions. This strengthens the notion that the uncertainties in these measurements are likely underestimated. Every calibrated GBS field includes a distribution of sources that appear slightly brighter and slightly fainter than their Transient Survey counterparts. The brightest Strong candidates all display a positive $\dot{f}/{f}_{\mathrm{trans}}$, while the fainter Strong candidates all show significant brightness decreases between the GBS and Transient Survey eras. This may indicate an intrinsic difference in the underlying accretion process for different types of protostars. This trend, however, is only based on small number statistics.

Three factors contribute to the sensitivity of $| \dot{f}/{f}_{\mathrm{trans}}|$: the uncertainties in the GBS and Transient Survey peak brightness measurements, the uncertainty in the calibration, and the time lag between the two data sets for each field (see ${\rm{\Delta }}t$ in Table 1). In Table 5, we present both the sensitivity limit and the typical sensitivity in $| \dot{f}/{f}_{\mathrm{trans}}|$ for each field. To calculate the sensitivity limit for each field, we set Equation (5) (replacing $\bar{R}$ with the mean FCF and ${\sigma }_{R}$ with ${\sigma }_{\mathrm{FCF}}$) equal to 6 (the threshold at which we define a source to be a variable candidate) and solve for the GBS peak brightness in terms of the Transient Survey peak brightness. Then, we substitute the result into Equation (6) for each source and find the minimum allowed $| \dot{f}/{f}_{\mathrm{trans}}|$ for a variable candidate in each field. For the typical $| \dot{f}/{f}_{\mathrm{trans}}|$ sensitivity, we report the median $| \dot{f}/{f}_{\mathrm{trans}}|$ value calculated in this fashion. As the Transient Survey continues and we are able to observe longer time separations, our sensitivity will improve.

Table 5.  Current $\dot{f}/{f}_{\mathrm{trans}}$ Sensitivity Limits

GBS Field	${\rm{\Delta }}f$ Sensitivity Limit (%yr−1)	Typical ${\rm{\Delta }}f$ Sensitivity (%yr−1)
IC348-E	2.2	2.9
L1688-1	3.3	5.0
L1688-2	2.9	4.8
NGC1333-N	1.4	3.6
OMC1 tile4	2.0	4.3
OrionBN_450_S	3.7	7.4
OrionBS_450_E	3.8	17.7a
OrionBS_450_W	2.4	4.7
SerpensMain1	1.1	5.0
SerpensS-NE	1.0	4.6
SerpensS-NW	1.5	3.8

Note.

^aThe large uncertainties in OrionBS_450_E are due to complication in background structure subtraction executed by the GaussClumps algorithm (see Appendix A). These are, however, bright and isolated sources in this field that are well recovered and better represented by the $| \dot{f}/{f}_{\mathrm{trans}}|$ sensitivity limit in the adjacent column.

Download table as:  ASCIITypeset image

The relative flux calibration agreement between GBS fields associated with the same Transient Survey field is robust. There are 13 sources in total that overlap two GBS fields and a Transient Survey field. There are no cases where a source is considered to be a variable candidate in one GBS field and not the other. In total, we find seven sources with significantly discrepant GBS to Transient Survey peak brightness ratios, indicating they are candidates for variability. These sources are distributed around $\dot{f}/{f}_{\mathrm{trans}}=0 \% \,{\mathrm{yr}}^{-1}$ (Figure 8). Though we do not necessarily expect linear changes over 2–4 year timescales, assuming linearity, the average $| \dot{f}/{f}_{\mathrm{trans}}|$ of the Strong variable candidates over several year timescales using co-added, well-calibrated maps is $| 4.0| \% \,{\mathrm{yr}}^{-1}$, with a standard deviation of $2.7 \% \,{\mathrm{yr}}^{-1}$ (see Table 2 and Figure 8). We are not sensitive to short timescale, small fluctuations that may average out the signal over time (see Section 5.2 for a further example of the limitations of time averaging), so we expect there to be more variable candidates in these fields that can be identified at higher time resolutions and sensitivities. An analysis of variability within the Transient Survey data alone is addressed by D. Johnstone et al. (2017, in preparation).

The observed difference in the submillimeter flux of an embedded protostar undergoing a change in accretion rate is determined by the heating or cooling of the dusty envelope. As long as the dust temperature remains above ∼25 K, the sub-mm response will be approximately linear to this change in temperature. When the dust temperature is lower, the submillimeter response becomes much stronger. In equilibrium, the dust temperature, ${T}_{\mathrm{dust}}$, in most of the envelope is expected to vary with accretion luminosity, ${L}_{\mathrm{acc}}$, as ${T}_{\mathrm{dust}}\propto {L}_{\mathrm{acc}}^{1/4}$ (Johnstone et al. 2013). As Johnstone et al. (2013) note, however, this relationship will break down in the very outer envelope where significant heating comes from the external radiation field and therefore the temperature of the dust remains relatively fixed.

For small fractional changes in temperature, the accretion luminosity is approximated by

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Delta }}{L}_{{\rm{p}}}}{{L}_{{\rm{p}}}}\propto 4\times \displaystyle \frac{{\rm{\Delta }}{T}_{\mathrm{dust}}}{T}.\end{eqnarray} \tag{ 7 }$$

Assuming the submillimeter response is roughly proportional to the dust temperature, which itself is determined primarily by the accretion luminosity, a 4% change in the observed submillimeter flux corresponds to a ∼16% change in the accretion luminosity (accretion rate) of the central protostar. If the envelope temperature drops below ∼25 K, then the submillimeter response will be closer to one-to-one (i.e., a 4% change in the observed submillimeter flux corresponds to ∼4% change in accretion luminosity). Much work is still required to fully understand the relationship between the observed submillimeter flux and the accretion luminosity (see, e.g., Johnstone et al. 2013).

We find 7% of the known protostars in our observed fields display a typical 16% accretion variability over ∼3 years (using Equation (7)). In Figure 3 of Herczeg et al. (2017), the authors summarize expectation values of accretion variability over specific timescales, taking into consideration a model in which the variability is driven by large-scale gravitational instabilities in the disk (Vorobyov & Basu 2010) and a model in which smaller-scale magnetorotational instabilities are included Bae et al. (2014). The models of Vorobyov and Basu (2010) predict that ∼7% of protostars will undergo a 7%–8% change in accretion luminosity over 3 years of observations. Bae et al. (2014), however, predict that $\sim 7 \%$ of protostars will undergo a 40% change in accretion luminosity over this same time frame. Under the simple assumption that submillimeter brightness varies linearly as dust temperature, Equation (7) predicts our results to lie between the two models. We note, however, that this result is tempered by uncertainties in the relationship between changes in the submillimeter flux due to the protostellar luminosity, the reliability expected from the models of protostellar episodic accretion over few year timescales (as a detailed investigation relies on accurately tracing the physics of the inner disk), and the sensitivity of our source detection methods.

Future surveys studying accretion variability onto deeply embedded protostars would benefit from working in the far-infrared, near the peak of the protostellar envelope SED, in order to have a more linear relationship between observed brightness changes and the underlying accretion luminosity variation (see Johnstone et al. 2013). These surveys are likely to be undertaken by space observatories such as the Space Infrared Telescope for Cosmology and Astrophysics (Roelfsemaet al., submitted) or the Origins Space Telescope (Meixneret al. 2016), which will further benefit from the stability of space-based observations, provided that the calibration and dynamic range of the instruments is excellent. As these missions are still a decade or more away, an obvious first undertaking will be to use Herschel Space Observatory observations of nearby star-forming regions as a previous epoch, yielding a multi-decade delta in time. As this paper shows, however, there are significant challenges in collating two disparate data sets, and thus care will need to be taken in order to reach relative uncertainties at the 1%–2% level.

5.1. Previously Known Signatures of Variability

5.2. The Period of the Submillimeter Variable EC 53

In this paper, we investigate signatures of 850 μm submillimeter variability by comparing archival GBS fields, observed between 2012 and 2014, to JCMT Transient Survey fields, observed between 2015 December and 2017 March 1. We follow the data reduction and calibration procedures presented in Mairs et al. (2017) to self-consistently align and calibrate each set of observations from the two surveys individually, before bringing the co-added GBS images into relative flux calibration with the co-added Transient Survey images (see Section 3 and Figures 1 through 3). Using the source extraction algorithm Gaussclumps (Stutzki & Guesten 1990), we identify 175 independent bright (>200 mJy beam⁻¹), compact (effective radius <10'') emission objects that were well fit with Gaussian profiles, correlate them with known protostars (Class 0+I, Flat spectrum sources; see Section 2.3) and disks (Class II sources), and identify objects that are significant outliers in their GBS to Transient Survey average peak brightness ratio with respect to the other sources in the field (see Figures 4 through 6, Figure 7, and Equation (5)). Based on a visual inspection for compact structure present around the area of each outlying source in calibrated difference maps (the Transient Survey co-add subtracted from the GBS co-add), we define two categories for the quality of each variable candidate.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Our main results are summarized as follows:

• 1.  
  We have developed methods to robustly analyze repeated observations of an area of the sky for signatures of submillimeter variability. We identify 11 archival GBS fields that can be self-consistently flux calibrated and have significant overlap with the 8 Transient Survey fields (see Appendix A).
• 2.  
  Out of 175 independent compact, 850 μm emission sources, we find a total of 7 variable candidates, 4 of which are associated with known protostars.
• 3.  
  We classify 5 of the variable candidates as Strong and 2 as Extended (see Section 4), and we highlight one additional source, HOPS 383, as a Possible detection of variability. Two Strong sources (IRAS4A, IRAS 18270-0153), along with HOPS 383, have previously noted signatures of variability. There is one Strong variable candidate without a known protostellar or disk association (Bolo 40/PER-10). This is a good target for follow-up studies to identify signatures of very faint, deeply embedded protostars (see Table 2).
• 4.  
  The average flux change for the strong variable candidates is $| 4.0| \% \,{\mathrm{yr}}^{-1}$ with a standard deviation of $2.7 \% \,{\mathrm{yr}}^{-1}$ over 2–4 year timescales. Assuming the heating of the envelope is responsible for the changing luminosity, this corresponds to a change in accretion rate of $\sim 16 \%$ (see Section 5). The observed changes in flux for all sources are distributed around $0 \% \,{\mathrm{yr}}^{-1}$.
• 5.  
  Using the archival GBS data, we strengthen the detection of the submillimeter variable source EC 53 by adding a critical data point to its periodic light curve and determine its 850 μm period to be 567 ± 32 days. This value is consistent with the 543 day period previously reported in the K-band (see Figure 9, Yoo et al. 2017, and Hodapp et al. 2012).

Zoom In Zoom Out Reset image size

Throughout this work, we have developed methods to robustly analyze repeated observations of an area of the sky for signatures of submillimeter variability. The JCMT Transient Survey will continue through at least 2019 January, and as more data are collected, we will have the opportunity to continue this investigation to fainter sources and smaller levels of variability over longer timescales. Future directions include periodically co-adding sets of Transient Survey observations to construct very precise light curves (sacrificing temporal resolution for higher sensitivity), improving the flux calibration procedures such that more archival data can be included, comparing these observations with the simultaneous 450 μm data sets once the latter can be robustly reduced and calibrated, and applying these observational constraints to current simulations of variable protostellar accretion.

Steve Mairs was partially supported by the Natural Sciences and Engineering Research Council (NSERC) of Canada graduate scholarship program. Doug Johnstone is supported by the National Research Council of Canada and by an NSERC Discovery Grant. Gregory Herczeg is supported by general grant 11473005 awarded by the National Science Foundation of China. Andy Pon received partial salary support from a Canadian Institute for Theoretical Astrophysics (CITA) National Fellowship. Miju Kang was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF), funded by the Ministry of Science, ICT, & Future Planning (No. NRF-2015R1C1A1A01052160). J.-E. Lee was supported by the Basic Science Research Program through the National Research Foundation of Korea (grant No. NRF-2015R1A2A2A01004769) and the Korea Astronomy and Space Science Institute under the R&D program (Project No. 2015-1-320-18) supervised by the Ministry of Science and ICT.

The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain. The James Clerk Maxwell Telescope is operated by the East Asian Observatory on behalf of the National Astronomical Observatory of Japan, Academia Sinica Institute of Astronomy and Astrophysics, the Korea Astronomy and Space Science Institute, the National Astronomical Observatories of China, and the Chinese Academy of Sciences (Grant No. XDB09000000), with additional funding support from the Science and Technology Facilities Council of the United Kingdom and participating universities in the United Kingdom and Canada. The James Clerk Maxwell Telescope has historically been operated by the Joint Astronomy Centre on behalf of the Science and Technology Facilities Council of the United Kingdom, the National Research Council of Canada, and the Netherlands Organisation for Scientific Research. Additional funds for the construction of SCUBA-2 were provided by the Canada Foundation for Innovation. The identification number for the JCMT Transient Survey data used in this paper is M16AL001. The identification numbers for the archival Gould Belt Survey data used in this paper are MJLSG31, MJLSG32, MJLSG33, MJLSG38, and MJLSG41. The authors thank the JCMT staff for their support of the data collection and reduction efforts. This research has made use of NASA's Astrophysics Data System and the facilities of the Canadian Astronomy Data Centre operated by the National Research Council of Canada, with the support of the Canadian Space Agency. The authors would especially like to thank Chang Won Lee and Harriet Parsons for their useful insights and suggestions, along with the extended JCMT Transient Team²¹ for their support. This research used the services of the Canadian Advanced Network for Astronomy Research (CANFAR), which in turn is supported by CANARIE, Compute Canada, University of Victoria, the National Research Council of Canada, and the Canadian Space Agency. This research made use of APLpy, an open-source plotting package for Python hosted at http://aplpy.github.com, and matplotlib, a 2D plotting library for Python (Hunter 2007).

Facility: JCMT (SCUBA-2) (Holland et al. 2013).

Software: Starlink (Currie et al. 2014), Astropy (Astropy Collaboration et al. 2013), Python version 3.5, APLpy (Robitaille & Bressert 2012), Matplotlib (Hunter 2007).

For maximum consistency in relative flux calibration across all the regions, we select individual GBS fields and compare each one with its respective Transient Survey field independently (see Figures 10–17). For a given GBS field to be useful for this work, it must (1) have enough bright, compact sources to undergo self-calibration; and (2) have significant overlap with a Transient field such that they share at least three sources of interest and a flux calibration between the two data sets can be performed. The first criterion generally causes GBS fields with some Transient field overlap to fail. In order for a specific GBS field to undergo relative flux calibration, we require at least two compact sources brighter than 0.5 Jy beam⁻¹ that are observed to have consistent peak brightnesses with respect to one another across all observations. If two potential calibrator sources are noted but display discordant flux calibration factors, it is unclear which, if either, of the sources represents the more correct value without further investigation. Where we lack robust sources, we discard the field. Similarly, if only one potential calibrator source is noted, it is unclear whether or not that source is intrinsically variable and, again, we discard the field. A self-calibration failure is generally encountered when the field is very sparse, or when a field is inundated with complex, extended emission without isolated, well-defined point-like sources. In the latter case, the Gaussclumps algorithm has difficulties in defining source boundaries and properly subtracting background emission (see, e.g., the "Integral Shaped Filament" in the OMC 2–3 region in Figure 12; two GBS fields with the southern extension of the filament present along the noisy edge of their field of view were discarded).

Note that a clear example of Gaussclumps having difficulty subtracting the background is in the ORIONBS_450_E field. The error bars for many sources in Figure 2 (bottom left) are large, due to the inconsistency of the identified source boundaries. To appear in Figures 1 through 3, a source must be consistently observed in at least 2 GBS observations and at least 2 Transient Survey observations. Recall that we apply minimum peak brightness and maximum radius thresholds. Sources that have sizes that are near this threshold may be culled from many observations, yet detected in two or three, as the extraction algorithm fits a slightly different Gaussian model to the same region on different dates, attempting to separate clustered structure. The lack of robust observations increases the uncertainty of the peak measurement, but this is mitigated by well-defined, well-recovered sources in conjunction with the calculation of a weighted mean to derive the relative FCF.

Source recovery in particularly confused and blended emission regions is also complicated by the fact that the recovery of extended background emission can vary from observation to observation, especially near the edges of the map. Reducing SCUBA-2 data is a complex process. Best practices for the GBS data reduction were developed by Mairs et al. (2015) and extended to the Transient Survey by Mairs et al. (2017). Adopting the methodology of Mairs et al. (2017) in this paper, we robustly recover compact structure in exchange for less information at extended scales. As discussed in Section 3.1, we perform a two-stage reduction where we employ an external mask in the final stage to help constrain makemap's solution. These external masks are created individually for each field in the GBS and Transient Surveys, due to the different central coordinates of the overlapping fields. Where the fields from the two surveys overlap (see Figures 10 through 17), the masked structures are nearly identical. Outside of these regions, however, the fields have their own structure that can slightly affect how astronomical signal across the rest of the recovered image grows in the final iteration of makemap. This effect is generally insignificant across the majority of the image but can cause slight differences near the edges of the fields where the noise levels are higher.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In many instances, the edges of GBS fields overlap with the center of their associated Transient Survey field and the small differences in extended structure recovery create flux pedestals and negative bowling, which add to the uncertainty of a measured source in those regions (see Mairs et al. 2015 for more information). Since the telescope has to slow down and speed up near the edges, the time domain filtering has a different effect on structure near the edge of the map with respect to that at the center of the map. We have addressed this by inspecting each of the variable candidate sources in Table 2 in difference maps we constructed (see Figures 18 and 19), looking for indications of larger-scale residual structures. We expect that compact, truly varying sources would show significant, point-like structure, even in the midst of these extended regions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figures 18 and 19, we present thumbnail images extracted from the constructed difference maps (the GBS co-added data subtracted from the Transient Survey co-added data) for all variable candidates in Table 2. All of the Strong variable candidates show significant, compact structure in their respective difference maps.