# Submillimeter Continuum Variability in Planck Galactic Cold Clumps

Geumsook Park, Kee-Tae Kim, Doug Johnstone, Sung-Ju Kang, Tie T. Liu, Steve Mairs, Minho Choi, Jeong-Eun Lee, Patricio Sanhueza, Mika Juvela, et al.

► To cite this version:

Geumsook Park, Kee-Tae Kim, Doug Johnstone, Sung-Ju Kang, Tie T. Liu, et al. Submillimeter Continuum Variability in Planck Galactic Cold Clumps. The Astrophysical Journal Supplement, 2019, 242 (2), pp.27. 10.3847/1538-4365/ab1eae. hal-02397468

HAL Id: hal-02397468

https://hal.science/hal-02397468v1

Submitted on 26 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Submillimeter Continuum Variability in Planck Galactic Cold Clumps

Geumsook Park<sup>1</sup>, Kee-Tae Kim<sup>1</sup>, Doug Johnstone<sup>2,3</sup>, Sung-ju Kang<sup>1</sup>, Tie Liu<sup>1,4</sup>, Steve Mairs<sup>4</sup>, Minho Choi<sup>1</sup>, Jeong-Eun Lee<sup>5</sup>, Patricio Sanhueza<sup>6</sup>, Mika Juvela<sup>7</sup>, Miju Kang<sup>1</sup>, David Eden<sup>8</sup>, Archana Soam<sup>9</sup>, Julien Montillaud<sup>10</sup>, Gary A. Fuller<sup>11</sup>, Patrick M. Koch<sup>12</sup>, Chang Won Lee<sup>1,13</sup>, Dimitris Stamatellos<sup>14</sup>, Jonathan Rawlings<sup>15</sup>, Gwanjeong Kim<sup>16</sup>, Chuan-Peng Zhang<sup>17</sup>, Woojin Kwon<sup>1,13</sup>, and Hyunju Yoo<sup>1</sup>

*Received 2018 December 16; revised 2019 April 29; accepted 2019 April 30; published 2019 June 17*

## Abstract

In the early stages of star formation, a protostar is deeply embedded in an optically thick envelope such that it is not directly observable. Variations in the protostellar accretion rate, however, will cause luminosity changes that are reprocessed by the surrounding envelope and are observable at submillimeter wavelengths. We searched for submillimeter flux variability toward 12 Planck Galactic Cold Clumps detected by the James Clerk Maxwell Telescope (JCMT)-SCUBA-2 Continuum Observations of Pre-protostellar Evolution (SCOPE) survey. These observations were conducted at 850 µm using the JCMT/SCUBA-2. Each field was observed three times over about 14 months between 2016 April and 2017 June. We applied a relative flux calibration and achieved a calibration uncertainty of ~3.6% on average. We identified 136 clumps across 12 fields and detected four sources with flux variations of ~30%. For three of these sources, the variations appear to be primarily due to large-scale contamination, leaving one plausible candidate. The flux change of the candidate may be associated with low- or intermediate-mass star formation assuming a distance of 1.5 kpc, although we cannot completely rule out the possibility that it is a random deviation. Further studies with dedicated monitoring would provide a better understanding of the detailed relationship between submillimeter flux and accretion rate variabilities while enhancing the search for variability in star-forming clumps farther away than the Gould Belt.

**Key words:** stars: formation – submillimeter: general – submillimeter: ISM – surveys

## 1. Introduction

A protostar gains mass by accreting material through a protostellar disk embedded in a circumstellar envelope (for a review, see e.g., Hartmann et al. 2016). Understanding the mass accretion of protostars is an essential component to characterizing their overall formation and evolution. The earliest formulations of star formation theory assumed a steady-state accretion model where the amount of mass gained by the protostar was constant over time (Shu 1977; Terebey et al. 1984). Many subsequent observational studies (e.g., Kenyon et al. 1990; Evans et al. 2009), however, indicate that protostellar luminosities are lower than predicted by these conventional models. One solution to this so-called “luminosity problem” is a variable protostellar accretion rate (often called episodic accretion), where bright outbursts occur over short timescales and the forming star spends most of its time in a

“quiescent” phase (Dunham et al. 2010; Dunham & Vorobyov 2012). The evolutionary lifetime of protostars, however, is still uncertain and refining the current estimates may also contribute to correcting this apparent discrepancy between the protostellar luminosity predicted by models and the current observations (e.g., Evans et al. 2009; McKee & Offner 2011). Ultimately, studies of the variability of accretion rates are critical in order to understand the physics of the circumstellar disk and how the mass is transported onto the protostar, itself.

In this study, we focus on detecting signs of episodic accretion in the earliest stages of star formation. The majority of accretion variability observations have so far been carried out in the evolved stages of pre-main-sequence stars (e.g., Kóspál et al. 2007; Aspin et al. 2009; Caratti o Garatti et al. 2011; Covey et al. 2011; Fischer et al. 2012; Reipurth et al. 2012). EX Lupi (e.g., Herbig 2008; Aspin et al. 2010) and FU Orionis (e.g., Herbig 1977; Hartmann & Kenyon 1996) sources

<footnote n="1"><sup>1</sup> Korea Astronomy and Space Science Institute, 776 Daedokdae-ro, Yuseong-gu, Daejeon 34055, Republic of Korea; pgs@kasi.re.kr</footnote>

<footnote n="2"><sup>2</sup> NRC Herzberg Astronomy and Astrophysics, 5071 West Saanich Rd., Victoria, BC, V9E 2E7, Canada</footnote>

<footnote n="3"><sup>3</sup> Department of Physics and Astronomy, University of Victoria, Victoria, BC, V8P 5C2, Canada</footnote>

<footnote n="4"><sup>4</sup> East Asian Observatory, 660 North A'ohoku Pl., University Park, Hilo, Hawaii 96720, USA</footnote>

<footnote n="5"><sup>5</sup> School of Space Research, Kyung Hee University, 1732, Deogyong-Daero, Giheung-gu, Yongin-shi, Gyunggi-do 17104, Republic of Korea</footnote>

<footnote n="6"><sup>6</sup> National Astronomical Observatory of Japan, National Institutes of Natural Sciences, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan</footnote>

<footnote n="7"><sup>7</sup> Department of Physics, P.O. Box 64, FI-00014, University of Helsinki, Finland</footnote>

<footnote n="8"><sup>8</sup> Astrophysics Research Institute, Liverpool John Moores University, IC2, Liverpool Science Park, 146 Brownlow Hill, Liverpool L3 5RF, UK</footnote>

<footnote n="9"><sup>9</sup> SOFIA Science Centre, USRA, NASA Ames Research Center, MS-12, N232, Moffett Field, CA 94035, USA</footnote>

<footnote n="10"><sup>10</sup> Institut UTINAM, CNRS UMR 6213, OSU THETA, Université Bourgogne Franche-Comté, France</footnote>

<footnote n="11"><sup>11</sup> Jodrell Bank Centre for Astrophysics, School of Physics and Astronomy, The University of Manchester, Oxford Rd., Manchester M13 9PL, UK</footnote>

<footnote n="12"><sup>12</sup> Institute of Astronomy and Astrophysics, Academia Sinica. 11F of Astronomy-Mathematics Building, AS/NTU No. 1, Section 4, Roosevelt Rd., Taipei 10617, Taiwan</footnote>

<footnote n="13"><sup>13</sup> University of Science and Technology, Korea (UST), 217 Gajeong-ro, Yuseong-gu, Daejeon 34113, Republic of Korea</footnote>

<footnote n="14"><sup>14</sup> Jeremiah Horrocks Institute for Mathematics, Physics & Astronomy, University of Central Lancashire, Preston PR1 2HE, UK</footnote>

<footnote n="15"><sup>15</sup> Department of Physics and Astronomy, University College London, Gower St., London, WC1E 6BT, UK</footnote>

<footnote n="16"><sup>16</sup> Nobeyama Radio Observatory, National Astronomical Observatory of Japan, National Institutes of Natural Sciences, Nobeyama, Minamimaki, Minamisaku, Nagano 384-1305, Japan</footnote>

<footnote n="17"><sup>17</sup> National Astronomical Observatories, Chinese Academy of Sciences, Beijing, 100012, People's Republic of China</footnote>
could be classical examples of episodic accretion occurring after the deeply embedded phase. The spectacular observational change in optical brightness for these objects is about a factor of 10 or more and lasts for several months to decades. Recently, however, a few outbursts from deeply embedded protostellar objects have been reported (e.g., Safron et al. 2015; Hunter et al. 2017; Yoo et al. 2017).

The mass accretion rates are high in the early stages of protostellar evolution (e.g., Whitworth & Ward-Thompson 2001; Schmeja & Klessen 2004), so we expect the accretion variability to be more significant than during later stages. Direct observations at optical or near-infrared (near-IR) variability of protostellar systems (star(s)/disk(s)) are, however, very challenging because these systems are heavily embedded in optically thick, dense envelopes. Thus, indirect observations at submillimeter wavelengths are now being explored. Johnstone et al. (2013) analyzed the flux variability of a protostellar envelope caused by outbursts of the central source using far-IR and submillimeter continuum emission. The model suggests that mid- to far-IR observations would be ideal to detect variability changes over timescales of hours to days. The study also revealed that detecting variability at submillimeter wavelengths should be achievable, although variations occur over longer timescales of approximately one month. There are, indeed, some examples of known flux variations that were detected in submillimeter continuum emission: about a factor of 2 flux increase at 350 µm and 450 µm toward the Class 0 source HOPS 383 (Safron et al. 2015). More recently, Mairs et al. (2018) reported that HOPS 358 now shows a strong, declining light curve over the course of 16 months. The Class I protostar, EC 53, displayed a 50% flux increase at 850 µm (Yoo et al. 2017). Furthermore, Hunter et al. (2017) reported high-mass protostellar system NGC 6334I-MM1 with a factor of 4.2 increase in 870 µm continuum interferometric flux and a 30% increase in the submillimeter single-dish flux.

As one of the large programs at the East Asian Observatory's James Clerk Maxwell Telescope (JCMT), the JCMT Transient Survey (Herczeg et al. 2017) has been designed to search for this type of long-term variability in submillimeter dust emission surrounding deeply embedded protostars in eight nearby star-forming regions within ~500 pc from the Sun. The Transient Survey is the only monitoring survey performed at submillimeter wavelengths. The first major results were released after 1.5 yr, and the team has also used archival data to identify variability over a timescale of ~5 yr (Johnstone et al. 2018; Mairs et al. 2017a). They found that ~10% of deeply embedded protostars display varying flux at the level of 5%–10% per year. However, these nearby regions are mostly forming low-mass stars.

The SCUBA-2 Continuum Observations of Pre-protostellar Evolution (SCOPE) survey with the JCMT (Liu et al. 2018a) has observed ~1200 Planck Galactic cold clumps (PGCCs; Planck Collaboration et al. 2016), which were selected in wide ranges of Galactic longitudes and latitudes. The SCOPE sample was biased to high column density PGCCs of $N_{H_2} > 1 \times 10^{21} \text{cm}^{-2}$ (in Planck measurements), but also included randomly selected lower column density PGCCs ($>5 \times 10^{20} \text{cm}^{-2}$ in Planck measurements). For about 3/5 of the SCOPE sample, physical properties are given in Planck Collaboration et al. (2016): (1) about 70% among them are concentrated within 1 kpc while

the others are widely distributed at up to ~8 kpc, with an average angular size of ~8'; (2) the mass range is from 0.1 $M_\odot$ to $10^5 M_\odot$ (see Figure 2 of Liu et al. 2018a and Figure 1 of Eden et al. 2019 for detailed distributions). Therefore, the SCOPE sample contains diverse clumps in different Galactic environments, from low-mass to high-mass star-forming regions at various distances from the Sun. In addition, to obtain deep images of high-mass star-forming regions as well as to detect large flux variation events, the SCOPE survey observed some (<30) PGCCs, which are composed of multiple substructures, on three separate occasions. We note that the SCOPE survey looks at more distant clumps than the JCMT Transient Survey and, thus, they are more likely to contain groups of protostars rather than individuals. When accretion variability is detected, therefore, the flux can be diluted if the event originates in a single protostar as the beam contains many protostars. Nevertheless, as shown in previous studies by Mairs et al. (2017b) and Johnstone et al. (2018), we expect to uncover flux variability at about the 10%-level or larger.

In this paper, we examine the flux variability from 12 PGCC fields in the first quadrant of the Galactic plane using the SCOPE survey data that are described in Section 2. The data reduction, including calibration and clump identification, is presented in Section 3. We follow the procedures of the Transient Survey team (e.g., Mairs et al. 2017b; Johnstone et al. 2018) with appropriate modifications. We present the results of the examination of flux variability toward identified clumps in Section 4 and discuss possible candidates of flux variation in Section 5. Section 6 summarizes the main results.

## 2. Data

The SCOPE survey mapped approximately 1200 PGCCs (Planck Collaboration et al. 2016) at 850 µm using SCUBA-2, the submillimeter continuum imaging instrument (Holland et al. 2013) at the 15 m JCMT. The survey was begun in 2015 December and was completed in 2017 July. Each map is about 12' in diameter, and the main beam size (FWHM) of JCMT/SCUBA-2 is 14"/1 at 850 µm (Dempsey et al. 2013). Each field was observed under grade 3/4 weather conditions with zenith opacities at 225 GHz between 0.1 and 0.15. The first released data were obtained by filtering out scales larger than 200" in order to remove the effects of the atmosphere, which is bright in the submillimeter regime. The pixel size is 4". The applied flux conversion factor (FCF) is 554 Jy $\text{pW}^{-1}\text{beam}^{-1}$ (Liu et al. 2018a), which is slightly higher than the usual FCF of 537 Jy $\text{pW}^{-1}\text{beam}^{-1}$ given by Dempsey et al. (2013). Considering this research is a part of the SCOPE survey, we initially adopted the FCF values calculated from Liu et al. (2018a) in order to keep the consistency of the data. Nevertheless, since we performed a relative flux calibration (see Section 3 for details), the absolute flux calibration is not essential for this study. This survey provides ~20 times higher angular resolution images compared to the Planck 353 GHz (~850 µm) data (~5'). Thus, complex substructures in the PGCCs can be resolved in the JCMT images which could not be resolved by Planck (e.g., see Figure 1). Liu et al. (2018a) presents a detailed description of the survey, and Eden et al. (2019) provides information of the first data release and the catalog of compact sources resolved with the JCMT.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Co-added images for the 12 fields observed by the SCOPE survey. Each image is cropped using a circle with a radius of 370". The field name is displayed at the top left of each image. The color bar is shown in a linear scale, ranging from -30 mJy beam⁻¹ (black) to the value in parentheses at the top right of each image. White circles are marked to help to locate the outliers described in Section 4, and the number assigned to each circle is from Table 3.</figcaption></figure>

In this study, we selected 12 PGCC fields in the first quadrant of the Galactic plane that are moderately bright and contain a relatively large number of clumps.<fnref n="18" /> PGCCs are written using the acronym “PGCCs” in the text. These regions span the Galactic longitude range of $14° < l < 36°$ and are located at heliocentric distances from ~1.5 to 17 kpc (Table 1). The three observations of each field were not carried out with a regular cadence and, therefore, had intervals spanning three weeks to 13 months. The total exposure time to complete each epoch is 15.4 minutes on average, and the median and maximum of exposure times per pixel are ~55 and ~200 s, respectively. Each image was smoothed with a Gaussian kernel of 8" FWHM (twice the pixel size) to reduce pixel-to-pixel noise. Thus, the final images shown in this paper have an angular resolution of 16"/2 FWHM after smoothing.

### 3. Data Reduction

The default 850 µm absolute flux calibration produced by the data reduction pipeline at the JCMT yields a 5%–10% uncertainty in pointlike calibrator sources over weather bands 1 through 4 (Dempsey et al. 2013; Mairs et al. 2017b). Therefore, to detect a 3σ_rms change in the peak flux of a source, the brightness variation would need to be at least 15%–30%. Simulations (e.g., Bae et al. 2014; Vorobyov & Basu 2015) as well as JCMT Transient Survey observations (Mairs et al. 2017a; Johnstone et al. 2018), however, suggest that less dramatic flux variations are more common. In order to increase detection reliability, it is advantageous to calibrate the flux in a relative sense using the method presented by Mairs et al. (2017b). In this way, it is possible to reduce the (relative) flux uncertainty to 2%–3%, which allows for statistically significant measurements of ~6%–10% flux changes.

In our implementation of the relative flux calibration scheme, we restricted the sources with high (>25) signal-to-noise ratios (S/Ns; see Section 4.1 for details). However, unlike the Transient Survey procedure, we did not require that the sources are compact. In comparing the source fluxes of different

<footnote n="18">¹⁸ The definition of “clump” is ambiguous. In this paper, SCOPE clumps resolved in the JCMT images are at various distances (see Table 1) and can contain substructures that are visible at higher resolution. The SCOPE clumps shown in this paper encompass masses from tens of solar masses to thousands of solar masses, spanning the range of cores to clouds. For simplicity, we refer to all these objects as clumps.</footnote>
<table><caption>**Table 1**
Fields and Epochs</caption><thead><tr><th rowspan="2">Field</th><th colspan="2">Central Position<sup>a</sup></th><th rowspan="2">Three Epochs<br>(yyyy mm dd)</th><th colspan="2">Time Intervals<sup>b</sup></th><th rowspan="2">Distance(s)<sup>c</sup><br>(kpc)</th></tr><tr><th>(h:m:s)</th><th>(d:m:s)</th><th>(day)</th><th>(day)</th></tr></thead><tbody><tr><td>G14.14−0.55</td><td>18:18:11.50</td><td>−16:55:29.05</td><td>2016 Apr 10</td><td>2017 May 10</td><td>2017 May 27</td><td>1.5</td></tr><tr><td>G14.47−0.20</td><td>18:17:31.80</td><td>−16:28:00.46</td><td>2016 Apr 9</td><td>2017 May 11</td><td>2017 Jun 2</td><td>3.1 (11.5)</td></tr><tr><td>G14.71−0.19</td><td>18:17:59.80</td><td>−16:14:41.16</td><td>2016 Apr 9</td><td>2017 May 10</td><td>2017 Jun 2</td><td>3.1</td></tr><tr><td>G15.61−0.48</td><td>18:20:48.40</td><td>−15:35:41.29</td><td>2016 Apr 10</td><td>2017 May 11</td><td>2017 Jun 2</td><td>1.8 and 16.9</td></tr><tr><td>G23.68 + 0.57</td><td>18:32:23.20</td><td>−07:57:39.50</td><td>2016 Apr 11</td><td>2017 May 10</td><td>2017 Jun 3</td><td>5.8</td></tr><tr><td>G23.97 + 0.51</td><td>18:33:09.20</td><td>−07:43:48.16</td><td>2016 Apr 11</td><td>2017 May 12</td><td>2017 Jun 4</td><td>5.8</td></tr><tr><td>G24.04 + 0.26</td><td>18:34:10.40</td><td>−07:47:05.86</td><td>2016 Apr 11</td><td>2017 May 10</td><td>2017 Jun 2</td><td>7.8</td></tr><tr><td>G24.49−0.52</td><td>18:37:48.10</td><td>−07:44:45.61</td><td>2016 Apr 11</td><td>2017 May 12</td><td>2017 Jun 2</td><td>11.3</td></tr><tr><td>G25.68−0.14</td><td>18:38:39.10</td><td>−06:30:49.20</td><td>2016 Apr 11</td><td>2017 May 9</td><td>2017 May 27</td><td>10.2 (7.4)</td></tr><tr><td>G26.17 + 0.13</td><td>18:38:34.70</td><td>−05:57:20.53</td><td>2016 Apr 11</td><td>2016 Aug 30</td><td>2017 Jun 4</td><td>7.6</td></tr><tr><td>G33.72−0.02</td><td>18:52:55.20</td><td>+00:41:26.00</td><td>2016 Apr 12</td><td>2016 Jul 22</td><td>2017 May 27</td><td>6.5 (2.2)</td></tr><tr><td>G35.49−0.31</td><td>18:57:12.90</td><td>+02:07:52.72</td><td>2016 Apr 13</td><td>2016 Jun 7</td><td>2017 May 27</td><td>2.7 (3.2 and 10.3)</td></tr></tbody></table>

**Notes.**

a Equatorial coordinates, R.A. and decl. (J2000).

b Time intervals between the first and second epochs and between the second and third epochs.

c Distances are obtained from (Urquhart et al. 2018, see also references therein). For fields having clumps at various distances, we give the distance of the majority of clumps along with the value(s) of the minority in parenthesis, or, if they are almost equal numbers, two values with the conjunction “and.”

epochs, we checked that the locations of most peaks remain within the nominal 2"–6" uncertainty of the JCMT pointing; a typical difference is ~4". Further, we applied an image registration technique. We used the IDL/SUBREG procedure<fnref n="19" /> and derived offsets ($\Delta_{ra}$ and $\Delta_{dec}$) between different images. The two-dimensional offsets ($\Delta = \sqrt{(\Delta_{ra})^2 + (\Delta_{dec})^2}$) of our SCOPE fields have a mean of 4" and a standard deviation of 2". This is consistent with our previous, manual inspection. For these reasons, in this research, the different peak positions in a given clump area within the beam size, were assumed to originate from the same source.

As with Mairs et al. (2017b), the peak flux values were compared in order to determine the relative calibration. For point sources, calibrating to the peak flux allows a single number (the relative flux calibration) to be determined for each epoch, independent of the many underlying physical aspects responsible for the original calibration uncertainty. Sources of calibration uncertainty include: a poor measurement of the sky opacity, changed throughput of the instrument, or a slight focus offset, the latter of which will contribute to a change in the observed beam shape. For extended sources, determining the relative calibration using only the peak flux introduces an additional level of uncertainty since changes to the underlying beam profile also produce changes in the expected flux of the source. Despite this complication, Mairs et al. (2017b) and Mairs et al. (2015) found that the peak flux of bright sources embedded in extended emission are well-recovered and consistent for data reduction methods similar to those used in this study. Furthermore, we derived a robust uncertainty associated with the relative flux calibration factor (RFCF; calculated below) as an additional check on the validity of the process.

In every observed field, each epoch was calibrated individually and co-added to produce a deep, averaged image (Figure 1). To achieve this, the PICARD package (Gibb et al. 2013) found in the STARLINK software (Currie et al. 2014) was used. Although each co-added image was made by combining three epoch images, the individual images were not very

accurately aligned. Therefore, we used the co-added image only for the clump identification without getting into the details of alignment. The peak flux density values used in the following analysis were obtained from the individual epochs.

We identified submillimeter clumps in the co-added images with the CLUMPFIND algorithm (Williams et al. 1994), provided by STARLINK’s CUPID package (Berry et al. 2007), considering an rms noise level described in Section 3.1. There are several parameters to be set, such as “FwhmBeam,” “MinPix,” “MaxBad,” and “Tlow.”<fnref n="20" /> During the implementation, resultant clumps containing fewer pixels than the area corresponding to the beam size (< MinPix) were discarded. In addition, we excluded any clump if its peak is located beyond 370" from the central position of each map.<fnref n="21" /> Information regarding the structure in each field, along with the derived relative calibration factors, are listed in Table 1.

The relative flux calibration using the SCOPE data started by assuming that none of the clumps are variable. We found stable calibrator sources by an iterative method. From the relative flux calibration derived using the stable sources, we achieved a sensitivity that is sufficient to robustly detect a 10% flux variation (see Section 3.2 for details). Then, we examined whether non-calibrator sources are outliers and tested their significance with respect to the observational uncertainty.

### 3.1. Measuring the Flux: Step 1

A robust rms noise measurement is important not only for identifying clumps but for assessing the significance of their flux variability. However, the rms noise level of a constant

<footnote n="20">²⁰ “FwhmBeam” defines the FWHM size of the JCMT beam in pixels, which corresponds to 4.05 for our final images. “MinPix” is the smallest number of pixels which a clump can have; we used a value of 13 as that corresponds to the area of a circle with a diameter equal to the (post-smoothing) beam FWHM. “MaxBad” is the maximum fraction of blank pixels that can be contained in a clump, which is set to zero. “Tlow” defines the lowest contour level to consider; we use 3×rms noise. A detailed description of the parameters is given at http://starlink.eao.hawaii.edu/docs/sun255.htx/sun255ss5.html.</footnote>

<footnote n="21">²¹ This corresponds to the radius of the images in Figure 1. The maps are shaped like an uneven circle of which the radius extends from ~6'5 to ~8'. Near the edges of the images, the fields were much less exposed and the coverage is uneven from epoch to epoch.</footnote>

<footnote n="19">¹⁹ http://www.stsci.edu/~mperrin/software/sources/subreg.pro</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Top: example of the rms noise as a function of exposure time. Black diamonds, green crosses, and red pluses represent (in order) the three epochs observed of the G14.14–0.55 field. For the rms noise calculation, we only used bins where there are more than 200 pixels (dashed gray line in the bottom panel). The best-fit noise profiles ~1/√t, where t is the exposure time, are marked. Their colors match the data points. Bottom: the number of pixels as a function of exposure time.</figcaption></figure>

velocity (CV) daisy observation,<fnref n="22" /> (the mode we employed in the SCOPE survey) is not uniform over the entire field. Since SCUBA-2 generates a map of the exposure time for each mapping field, we were able to use this map to characterize the rms noise at different positions in the field. We measured the rms noise level as a function of the exposure time in areas with no astronomical signal using data from each epoch. The rms noise levels showed gradual changes (almost flat) at exposure times larger than ~50 s and increase sharply at shorter exposure times (see Figure 2). By design, most pixels in the latter case are located near the edge of the uncropped images, so the data points with exposure times shorter than ~50 s are insignificant for our analysis. We generated a best-fit noise profile for each epoch (curves in the top panel of Figure 2) using a simplified equation of the expected noise level (y) y = c₁ + c₂/√t, where t indicates exposure time. In the exposure time range of 50–200 s, we took the average of the best-fit noise profiles of the individual epochs. The average noise level was then scaled down by a factor of √3 to account for the co-adding of the three epochs. Finally, this value was used to identify significant clumps in the co-added image. Though each image has the same exposure time, the data quality also depends on the amount of precipitable water vapor in the sky during the observations as well as on the elevation of the field. As shown in Figure 2, however, the data points over the three epoch are consistent, implying that the data quality is comparable from epoch to epoch. We measured the rms noise values for each of the 12 co-added images in order to perform clump identification. For the 12 co-added images, the averaged

mean value of the resultant rms noise levels for finding clumps is ~4 mJy beam⁻¹. In a single epoch image, the rms noise level reaches ~8 mJy beam⁻¹ in the central area with the longest exposure time.

## 3.2. Measuring the Flux: Step 2

We measured the peak flux $F_e(i)$ for each clump $i$ and epoch $e$. We denote the mean peak flux over the three epochs as $F_m(i)$. The peak flux measurements are robust as the 8" Gaussian smoothing mitigates pixel-to-pixel noise variations and the peak position uncertainty from epoch to epoch is less than beam size. In addition, we selected clumps with $F_m(i) \ge 250$ mJy beam⁻¹ for this analysis, which is ~25 S/N in a single epoch (noise ~10 mJy beam⁻¹). To find stable calibrator sources for relative flux calibration, we first assumed that all clumps are not variable. In each epoch $e$, we derived a RFCF as follows:

$$RFCF = \frac{\sum_{c=0}^{n_c} F_e(c)/F_m(c)}{n_c}, \quad (1)$$

where $c$ denotes a calibrator, and $n_c$ is the number of calibrators per field. Each epoch image was divided by its RFCF in order to calibrate the images relative to one another. From these relative flux calibrated images, we remeasured the peak fluxes in each epoch and compared the standard deviation, $\sigma_{std,meas}(i)$, of the clump fluxes with a fiducial standard deviation model, $\sigma_{std,fid}(i)$. The fiducial standard deviation model characterizes the uncertainty in a relative flux calibrated image based on the rms noise ($\sigma_{rms}(i)$) and the relative flux calibration uncertainty itself ($u_{cal}$; see Johnstone et al. 2018 for further details). $\sigma_{std,fid}(i)$ is calculated as follows:

$$\sigma_{\text{std,fid}}(i) = \sqrt{\sigma_{\text{rms}}(i)^2 + (u_{\text{cal}} \times F_{\text{m}}(i))^2}, \quad (2)$$

where $u_{cal}$ is

$$u_{\text{cal}} = \sqrt{\frac{\sum_{c=1}^{n_c} \sigma_{\text{std,meas}}(c)^2 / F_m(c)^2}{n_c - 1}} \quad (3)$$

Here, $\sigma_{rms}(i)$ is the mean value of the three epoch noise levels shown in Table 3, and $u_{cal}$ is given in the last column of Table 2.

The relative calibration steps were repeated using a clipping process to identify a set of stable calibrators. After applying the relative flux conversions for each epoch, we compared the expected uncertainty for each source ($\sigma_{std,fid}(i)$) with the measured value ($\sigma_{std,meas}(i)$). As discussed in Section 4.1, with only three measurements, we expected $\sigma_{std,meas} < 1.7 \times \sigma_{std,fid}$, which corresponds to a 95% of confidence level if there is no intrinsic variability. The numbers of identified clumps, calibrator sources, outliers, the RFCF at each epoch, and $u_{cal}$ are listed in Table 2 (see Section 4.1 for details on the outliers). Figure 3 shows histograms of the normalized RFCFs (normalized to the first epoch) and associated uncertainties. The normalized RFCFs were used to moderate the effects of small number statistics. The applied RFCFs were within the nominal flux calibration uncertainty of SCUBA-2 data at 850 µm (Dempsey et al. 2013; Mairs et al. 2017b). The median relative calibration uncertainty ($u_{cal}$) was found to be ~3.6%, which is slightly higher than what the Transient Survey team achieved (~2%). This slight increase in the relative calibration uncertainty is primarily due to two effects: the lower brightness

<footnote n="22">²² http://www.eaobservatory.org/jcmt/instrumentation/continuum/scuba-2/observing-modes/</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Histograms of the RFCF and the relative calibration uncertainty, $u_{cal}$, from Table 2. Values in the RFCF histogram were normalized to the first epoch for each field, so only RFCFs derived for the other epochs are counted.</figcaption></figure>

<table><caption>**Table 2**
Number of Clumps Found and Relative Calibration Information</caption><thead><tr><th rowspan="2">Field</th><th rowspan="2">All Clumps Found<br>>250 mJy beam<sup>-1</sup></th><th rowspan="2">Calibrators</th><th rowspan="2">Outliers</th><th colspan="3">RFCF at Each Epoch</th><th rowspan="2">u<sub>cal</sub><br>(%)</th></tr><tr><th>First</th><th>Second</th><th>Third</th></tr></thead><tbody><tr><td>G14.14-0.55</td><td>30</td><td>28</td><td>2</td><td>1.005</td><td>0.981</td><td>1.014</td><td>3.6</td></tr><tr><td>G14.47-0.20</td><td>19</td><td>19</td><td>0</td><td>1.014</td><td>0.953</td><td>1.033</td><td>4.5</td></tr><tr><td>G14.71-0.19</td><td>13</td><td>13</td><td>0</td><td>1.029</td><td>0.932</td><td>1.039</td><td>5.2</td></tr><tr><td>G15.61-0.48</td><td>6</td><td>6</td><td>0</td><td>0.994</td><td>0.995</td><td>1.010</td><td>1.8</td></tr><tr><td>G23.68+0.57</td><td>4</td><td>4</td><td>0</td><td>1.027</td><td>0.937</td><td>1.037</td><td>4.2</td></tr><tr><td>G23.97+0.51</td><td>3</td><td>3</td><td>0</td><td>1.005</td><td>0.988</td><td>1.010</td><td>2.8</td></tr><tr><td>G24.04+0.26</td><td>10</td><td>9</td><td>1</td><td>1.038</td><td>0.895</td><td>1.066</td><td>3.1</td></tr><tr><td>G24.49-0.52</td><td>4</td><td>4</td><td>0</td><td>1.025</td><td>0.981</td><td>0.995</td><td>4.9</td></tr><tr><td>G25.68-0.14</td><td>18</td><td>17</td><td>1</td><td>1.043</td><td>1.012</td><td>0.945</td><td>4.4</td></tr><tr><td>G26.17+0.13</td><td>6</td><td>6</td><td>0</td><td>1.085</td><td>0.902</td><td>1.013</td><td>3.4</td></tr><tr><td>G33.72-0.02</td><td>14</td><td>14</td><td>0</td><td>1.033</td><td>0.992</td><td>0.975</td><td>2.5</td></tr><tr><td>G35.49-0.31</td><td>9</td><td>9</td><td>0</td><td>1.056</td><td>0.948</td><td>0.996</td><td>1.7</td></tr></tbody></table>

limit used here for potential calibrators and the necessity to allow extended sources as calibrators. These differences from the Transient Survey are discussed in more detail below.

### 3.3. Differences in Methodology from the Transient Survey

We adopted the methods performed by the Transient Survey team to investigate peak flux changes over time. However, the SCOPE survey was not optimized for this type of work, so the following alterations to the Transient Survey methodology were applied.

First, our smoothing kernel size is slightly larger than that of the Transient Survey team (8" as opposed to 6"). Second, we used the CLUMPFIND algorithm while the Transient Survey team used GAUSSCLUMPS (Stutzki & Guesten 1990). Both of these algorithms provide almost the same results overall, but there are some differences in complex areas of a given map. Third, we applied a different set of criteria from the Transient Survey to select clumps from the catalogs obtained by using each algorithm. The Transient Survey team considered only sources which are very bright (>50 $\sigma_{rms}$) and compact (effective radius assuming a circular projected configuration <10"), and which appear in every epoch. Alternatively, we included less bright (>25 $\sigma_{rms}$) sources

and more extended sources. Fourth, the calibrator selection described above in this section differs from that of the Transient Survey team due to the difference in the number of bright sources. While we considered all the clumps to be potential calibrators at the beginning and then selected the invariable clumps, the Transient Survey team could be more selective as their fields contain many compact, bright clumps for the calibration such that the uncertainty from the noise was less than 5% (Mairs et al. 2017b). In spite of the differences in bright source selection for the relative flux calibration, the procedure presented in this study is sufficient to detect a flux variation of 10% ($3 \times u_{cal}$).

## 4. Results

### 4.1. Analysis of Peak Flux Measurement

We identified 136 clumps with $F_{m} \ge 250$ mJy beam$^{-1}$ across the 12 fields. Figure 4 shows the $\sigma_{std,meas}/\sigma_{std,fid}$ as a function of the mean peak flux density. Almost all clumps (132/136; marked with filled symbols in the figure) show little flux changes and are used as calibrators. Four outliers (open symbols) in three different SCOPE fields were detected.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** $\sigma_{\text{std,meas}}/\sigma_{\text{std,fid}}$ vs. $F_m$ for all identified clumps. Filled and open symbols are calibrators and outliers, respectively. A dashed line marks a threshold of 1.7 for calibrator sources (see the text for more details). The number assigned to each outlier is also marked.</figcaption></figure>

<table><caption>**Table 3**
Peak Flux of Outliers in 850 µm</caption><thead><tr><th rowspan="2">#</th><th rowspan="2">Field</th><th rowspan="2">Name<sup>a</sup></th><th colspan="2">Peak Position<sup>a</sup></th><th colspan="6">F<sub>c</sub> at Each Epoch<sup>b,c</sup></th><th rowspan="2">F<sub>m</sub><sup>c</sup></th><th rowspan="2">σ<sub>std,meas</sub><sup>c</sup></th><th rowspan="2">σ<sub>std,fid</sub><sup>c,d</sup></th><th rowspan="2">σ<sub>std,meas</sub><sup>c</sup> / σ<sub>std,fid</sub></th></tr><tr><th>R.A.(J2000)</th><th>Decl. (J2000)</th><th>First</th><th>Second</th><th>Third</th><th>F<sub>m</sub><sup>c</sup></th><th>σ<sub>std,meas</sub><sup>c</sup></th><th>σ<sub>std,fid</sub><sup>c,d</sup></th></tr></thead><tbody><tr><td rowspan="2">1</td><td>G14.14</td><td>G14.143</td><td>18:18:02.02</td><td>-16:53:57.09</td><td>317</td><td>(10)</td><td>397</td><td>(10)</td><td>378</td><td>(10)</td><td>364</td><td>42</td><td>16</td><td>2.6 (~0.1%)</td></tr><tr><td>-0.55</td><td>-0.508</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="2">2</td><td>G14.14</td><td>G14.210</td><td>18:18:29.89</td><td>-16:52:57.05</td><td>257</td><td>(11)</td><td>243</td><td>(10)</td><td>306</td><td>(10)</td><td>269</td><td>33</td><td>14</td><td>2.4 (~0.3%)</td></tr><tr><td>-0.55</td><td>-0.598</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="2">3</td><td>G24.04</td><td>G24.008</td><td>18:34:19.82</td><td>-07:50:29.89</td><td>281</td><td>(8)</td><td>213</td><td>(11)</td><td>293</td><td>(8)</td><td>263</td><td>43</td><td>12</td><td>3.6 (&lt;0.01%)</td></tr><tr><td>+0.26</td><td>+0.203</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="2">4</td><td>G25.68</td><td>G25.635</td><td>18:38:31.32</td><td>-06:32:53.20</td><td>297</td><td>(9)</td><td>230</td><td>(8)</td><td>251</td><td>(11)</td><td>259</td><td>34</td><td>15</td><td>2.3 (~0.5%)</td></tr><tr><td>-0.14</td><td>-0.126</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

**Notes.**

<sup>a</sup> Name contains each peak position in Galactic coordinates. It is determined from the epoch data with the highest peak flux.

<sup>b</sup> Values in parentheses are map noise levels.

<sup>c</sup> Units of mJy beam<sup>-1</sup>.

<sup>d</sup> See Equation (2). For each source, the mean of noise levels at three epochs and $u_{cal}$ listed in Table 2 are used.

<sup>e</sup> Values in parentheses indicate how reliable the explanation that the flux change is due to the observational uncertainty is (see Section 4).

Johnstone et al. (2018) searched for submillimeter variability in 1643 bright sources across eight star-forming regions using the first 18-month data of monthly observations obtained by the JCMT Transient Survey. Figure 2 of Johnstone et al. (2018) is similar to Figure 4 in this paper. Their results of $\sigma_{\text{std,meas}}/\sigma_{\text{std,fid}}$ are much more tightly constrained toward a value of 1. This is mainly due to their larger set of data (10–15 epochs) per region. EC 53, a known variable source in Serpens Main (Hodapp et al. 2012; Yoo et al. 2017), is an extreme outlier with a value of $\sigma_{\text{std,meas}}/\sigma_{\text{std,fid}} = 5.6$. We found no clump that shows similar, exceptional variability in our data.

To analyze how significant the outlier detections are, we constructed a simple statistical test of the null hypothesis that there is no variability beyond the flux changes due to the observational uncertainty. For 100,000 trials, we drew three peak values (to represent three epochs) at random from a normal distribution with a mean of a given peak value and a standard deviation of $\sigma_{\text{std,fid}}$. We measured $\sigma_{\text{std,meas}}$

from these three measurements, calculated $\sigma_{\text{std,meas}}/\sigma_{\text{std,fid}}$ for each trial, and examined the probability density function of $\sigma_{\text{std,meas}}/\sigma_{\text{std,fid}}$. We found that the probability density function depends only on the number of observational epochs. For the three epoch case, the $\sigma_{\text{std,meas}}/\sigma_{\text{std,fid}}$ distribution has a mean of 0.85 and a median of 0.83. $\sigma_{\text{std,meas}}/\sigma_{\text{std,fid}} \simeq 1.7$ and 3.1 give the cumulative probabilities of ~95% and $\gtrsim 99.99\%$, respectively. 100% minus the cumulative probability indicates the probability that the flux changes are simply due to the observational uncertainty. All four outliers in Figure 4 have $\sigma_{\text{std,meas}}/\sigma_{\text{std,fid}} \ge 2.3$, which corresponds to less than 0.5%. (This result is equivalent to identifying outliers at least $2.8\sigma$ from the mean in a normal distribution.) Therefore, they might be candidate variable sources.

The four outliers are listed in Table 3. The parameter of $\sigma_{\text{std,meas}}/\sigma_{\text{std,fid}}$ is a good, dimensionless indicator of flux variability. For the outliers, $\sigma_{\text{std,meas}}/\sigma_{\text{std,fid}}$ is between 2.3 and 3.6. Compared with EC 53, the outliers have much smaller
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Peak flux variations of the four selected outliers in Table 3. The number assigned to each outlier is written in the top right corner. Symbols and colors are described in Figure 4. The error bars represent $\sigma_{\text{std,fid}}$.</figcaption></figure>

values of $\sigma_{\text{std,meas}}/\sigma_{\text{std,fid}}$. In addition, all four outliers are relatively faint clumps (\~400 mJy beam$^{-1}$). The outliers are described below in more detail.

Figures 5 and 6 show the peak flux variations of the outliers at 850 µm, which are approximately six times the noise level. While it is difficult to define the variability timescale with a limited number of observations and an uneven observational cadence, we analyzed the trend of peak fluxes. Outliers 1 and 4 showed clear differences between the first and the two subsequent epochs. Outlier 2 showed no significant flux variations between the first two epochs separated by a year, but a sudden flux increase is detected between the second and third epochs separated by less than a month. Outlier 3 showed a clear difference in flux after the initial long time interval and also after the later, shorter time interval. The fluxes measured in the first and last epochs, however, were similar to one another. The JCMT Transient Survey found that the majority of variables uncovered have long-term (a number of years), rather than short-term, variations (monthly-to-yearly timescales; Johnstone et al. 2018), though only rare, extremely bright events allow the survey to uncover variations within individual epochs (Mairs et al. 2019). Further monitoring is required to confirm such short-term variations.

## 4.2. Large-scale Bias Check

Thus far, the technique we used in this paper is to compare the peak fluxes of different epochs for each clump after relative flux calibration. However, it is well known that submillimeter continuum map reconstruction often creates low-level, artificial, extended structures that may affect simple peak flux measurements. Such complications are more likely to arise across small crowded maps, such as those undertaken by SCOPE, as compared with the large, sparser Transient Survey fields. Thus, in this section we test whether the observed brightness variations from the four candidate variables are truly localized as expected for compact sources.

Thus, we aligned SCOPE images using the algorithm IDL/SUBREG mentioned in Section 3 and made difference maps using those epochs containing the minimum and maximum peak flux values. Figure 7 shows the flux difference maps of the four outlier candidates, zoomed in to localized areas of $2' \times 2'$. For each source, there are three panels: brightest and faintest epoch outlier images and their difference map. For Outlier 1, it appears that the majority of the flux change is located at the peak position. Therefore, we can confirm that the flux variation genuinely originates from the brightness of the localized source. On the other hand, for the other three sources (Outliers 2–4), between epochs the extended emission rises along with the peak flux increase. This can be seen most clearly in Outlier 4. For Outlier 2, there is a peaking-up trend above the background change by ~30 mJy beam$^{-1}$, which is only about half of the anticipated value from the peak flux analysis alone. For Outlier 3, there is an increase of about 60 mJy beam$^{-1}$ over the background change. However, this trend does not peak at the location of the source.

In summary, we find that three of the four candidate variables (Outliers 2–4) are closely associated with large-scale flux variations between epochs. As we do not expect to observe large variations in the brightness of an extended structure in star-forming regions and we are well aware of the likelihood of artificial large-scale structure created during the map-making process, we remove these three sources from any further analysis. Outlier 1 remains a “candidate” variable, although it is not particularly “robust” (see Section 4.1).
<figure><img src="image_8.png" /><figcaption>**Figure 6.** JCMT 850 µm images of the four outliers (top to bottom) listed in Table 3 for each of the three epochs (left to right). The images have been re-scaled by the relative calibration process described in Section 3. The outliers are marked using 15"-radius dotted circles in cyan. The assigned numbers are written in the right side of the circles in the leftmost images. The color scale corresponds to the 850 µm flux in mJy beam⁻¹.</figcaption></figure>

# 5. Discussion

## 5.1. *Variable Candidate Found in This Study*

Outlier 1 (G14.143–0.508) was found in the G14.14–0.55 field. The 28 invariable clumps in the field have

$\sigma_{\text{std,meas}}/\sigma_{\text{std,fid}} \approx 0.1–1.6$ with an average of 0.7, while the outlier has 2.6.

We investigated whether Outlier 1 shows signs of star formation, in which case the detected flux change could potentially be attributed to accretion variability. The clumps we
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Flux difference maps of the four outliers at JCMT 850 µm. For each outlier, there are three images: (left) the outlier-brightest epoch image, (center) the outlier-faintest epoch image, and (right) their difference map. White contours show a flux level of 230 mJy beam⁻¹ from the leftmost image. The outliers are marked using 15″-radius dashed circles in cyan.</figcaption></figure>

identify in this study were covered by the APEX Telescope Large Area Survey of the Galaxy (ATLASGAL; Schuller et al. 2009). We, therefore, searched for ATLASGAL clumps near the peak flux position of this outlier. Urquhart et al. (2018) derived the distances and physical properties (including evolutionary classification) of about 8000 ATLASGAL clumps in Galactic disk in the Galactic longitude from 5° to 60°. The ATLASGAL was conducted at 870 µm with a beam size of 19.2. The observing wavelength and beam size are comparable to ours. Note that the ATLASGAL survey has a typical noise level of 50–70 mJy beam⁻¹, which is one order of magnitude higher than that of the SCOPE survey. Outlier 1 is associated with ATLASGAL clump AGAL014.142–00.509 that has a

vLSR of 21.1 km s⁻¹. The kinematic distance was estimated to be 1.5 kpc (Urquhart et al. 2018).

This clump seems to be deeply embedded in an IR dark cloud filament. Urquhart et al. (2018) inferred Outlier 1 to be in a quiescent phase, because it is dark or weak at near- to far-IR wavelengths. The flux variation in a quiescent (seemingly starless) clump may sound contradictory. It can be explained, however, by the presence of at least one undetected heavily embedded (proto)star(s). For example, recent studies by Liu et al. (2018b) using single-dish telescopes and Contreras et al. (2018) using the Atacama Large Millimeter/submillimeter Array (ALMA) detected high accretion rates in massive quiescent cores, which are comparable to those found in
high-mass protostellar objects (see also Traficante et al. 2017). Also, the non-detection of an IR counterpart may be due to the sensitivity limits of existing mid/far-IR surveys (e.g., see Section 3.2.4 of Svoboda et al. 2016). Urquhart et al. (2019) found from a molecular line survey that ~35% of 29 quiescent ATLASGAL clumps in their sample have relatively high (30–50 K) rotation temperatures, suggesting the existence of internal heating protostellar object(s). Moreover, there are indeed several discoveries of compact bipolar molecular outflows in otherwise quiescent ATLASGAL clumps/cores (e.g., Feng et al. 2016; Tan et al. 2016; Pillai et al. 2019). The driving sources of the detected outflows were suggested to be massive protostars in the very early evolutionary stage. If Outlier 1 is a seemingly starless clump (in fact, not starless), high-resolution molecular observations could uncover that the clump is actually in the earliest stage of star formation. Therefore, a further investigation of this clump at higher angular resolution and sensitivity is required to uncover the embedded protostar(s).

The relationship between the bolometric luminosity and the envelope mass is useful for determining whether there is low- or high-mass star formation occurring (e.g., Molinari et al. 2008; Urquhart et al. 2014; Motte et al. 2018). Based on the luminosity (≈37 $L_\odot$) and the mass (≈23 $M_\odot$) of the associated clump, Outlier 1 is very likely related to low- or intermediate-mass star formation rather than high-mass star formation.

## 5.2. Comparison with Known Submillimeter Variable Sources

There are a few known submillimeter variable sources observed in low- and high-mass star-forming regions. They are relatively close (≤1.3 kpc), while Outlier 1 in this study appears to be slightly more distant. As an example of a low-mass, variable protostellar system, EC 53 is located at 436 pc (Ortiz-León et al. 2017), and it brightened from 960 to 1450 mJy beam⁻¹ at 850 µm in a 14′′6 single-dish beam (Yoo et al. 2017). If moved to a greater distance, the brightness of EC 53 will diminish significantly as the distance increases and, thus, it would fall below our sensitivity threshold. The area of an outburst associated with accretion variability is unresolved at the outliers' distances. Even if it were embedded in additional material, allowing the larger clump to be visible, the change in brightness of EC 53 would only be ~40 or 1 mJy at far distances of 1.5 or 10 kpc, respectively. These low flux variations are marginally detectable or undetectable levels in our observations. Alternatively, the massive protostellar system NGC 6334I-MM1 at 1.3 kpc (Chibueze et al. 2014; Reid et al. 2014) brightened by as much as 30% of the single-dish flux (≈0.3 × 65300 ≃ 19600 mJy beam⁻¹) at 850 µm in an 17′′5 beam (Sandell 1994; Hunter et al. 2017). In the same way, one expects to measure the increase of ~14700 and 300 mJy for an event like NGC 6334I-MM1 at distances of 1.5 and 10 kpc, respectively. These large variations would be easily detectable in our observations, but none of our SCOPE sources show such a dramatic change. The peak flux change of Outlier 1, thus, appears to be related to an event of an intermediate scale in terms of luminosity based on the EC 53 and NGC 3664I-MM1 case analyses.

This study suggests that long-term monitoring of distant star-forming regions with the JCMT is suitable for detecting submillimeter variability. If the variable candidate is confirmed, observations with higher angular resolution and sensitivity, using an interferometer such as the ALMA, will give us a better

understanding of their properties. Follow-up high-resolution observations will not only more easily detect any flux variability with little beam dilution but they will also reveal the embedded young stellar object(s) being responsible for variability events. As an example, the interferometric flux of NGC 6334I-MM1 increased by a factor of 4, which is much greater than a ~30% increase in the single-dish flux (Hunter et al. 2017).

## 6. Summary and Conclusion

We investigated the flux variations of submillimeter clumps in 12 PGCC fields in the first quadrant of the Galactic plane. The fields were observed three times over approximately 14 months using the JCMT/SCUBA-2, as part of the SCOPE survey. The survey was not optimized for detailed studies on flux variation and, therefore, the observations only cover three epochs with uneven time intervals. Nevertheless, taking into account the non-uniform noise distributions of the maps, we succeeded in examining relative flux changes by comparing the peak fluxes of identified clumps among epochs. We performed a relative flux calibration as described in Section 3, with a typical uncertainty of ~3.6%. In the 12 PGCC fields, we identified 136 clumps with mean peak flux densities larger than 250 mJy beam⁻¹ (≥25 S/N). From the peak flux analysis, we found four “outliers” that appear to vary in time. The average flux change at 850 µm is about 30%. We examined whether the peak flux changes of the outliers are well localized in the flux difference maps. Finally, only one (Outlier 1) of the four outliers is a plausible “candidate” and is not biased by the large scales. The detected flux variation in Outlier 1 may be related to episodic accretion events in the very early stage of low- or intermediate-mass star formation, considering a kinematic distance of 1.5 kpc, although we cannot completely exclude the possibility that it is a purely statistical random deviation. According to the existing observational data at near- to far-IR, the star-forming sign is less evident. However, the flux variability found here suggests an additional investigation of this region at higher angular resolution and sensitivity to uncover the deeply embedded protostar(s) in this clump. Further research employing long-term monitoring will be helpful not only to confirm our results but also to give a better understanding of the accretion processes in star formation.

We thank the anonymous referee for many valuable comments and suggestions which helped us improve the paper. D.J. is supported by NRC Canada and an NSERC Discovery Grant. T.L. is supported by the EACOA fellowship. J.E.L. is supported by the Basic Science Research Program through the National Research Foundation of Korea (grant No. NRF-2018R1A2B6003423) and the Korea Astronomy and Space Science Institute under the R&D program supervised by the Ministry of Science, ICT, and Future Planning. M.J. acknowledges the support of the Academy of Finland Grant No. 285769. M.K. is supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Science and ICT & Future Planning (No. NRF-2015R1C1A1A01052160). J.M. acknowledges the support by the Programme National “Physique et Chimie du Milieu Interstellaire” (PCMI) of CNRS/INSU with INC/INP co-funded by CEA and CNES. A.S. is supported by National Science Foundation Grant-1715876. C.W.L. is supported by the Basic Science Research Program through the National
Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (NRF-2019R1A2C1010851). W.K. was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF-2016R1C1B2013642). 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 Operation, Maintenance and Upgrading Fund for Astronomical Telescopes and Facility Instruments, budgeted from the Ministry of Finance (MOF) of China and administrated by the Chinese Academy of Sciences (CAS), as well as the National Key R&D Program of China (No. 2017YFA0402700). Additional funding support is provided by the Science and Technology Facilities Council of the United Kingdom and participating universities in the United Kingdom and Canada. The identification number for the JCMT SCOPE Survey under which the SCUBA-2 data were used in this paper is M16AL003.

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

Software: STARLINK(Currie et al. 2014), IDL.

## ORCID iDs

Geumsook Park https://orcid.org/0000-0001-8467-3736

Kee-Tae Kim https://orcid.org/0000-0003-2412-7092

Doug Johnstone https://orcid.org/0000-0002-6773-459X

Sung-ju Kang https://orcid.org/0000-0002-5004-7216

Tie Liu https://orcid.org/0000-0002-5286-2564

Steve Mairs https://orcid.org/0000-0002-6956-0730

Jeong-Eun Lee https://orcid.org/0000-0003-3119-2087

Patricio Sanhueza https://orcid.org/0000-0002-7125-7685

Mika Juvela https://orcid.org/0000-0002-5809-4834

Miju Kang https://orcid.org/0000-0002-5016-050X

Archana Soam https://orcid.org/0000-0002-6386-2906

Gary A. Fuller https://orcid.org/0000-0001-8509-1818

Patrick M. Koch https://orcid.org/0000-0003-2777-5861

Chang Won Lee https://orcid.org/0000-0002-3179-6334

Dimitris Stamatellos https://orcid.org/0000-0002-4502-8344

Gwanjeong Kim https://orcid.org/0000-0003-2011-8172

Chuan-Peng Zhang https://orcid.org/0000-0002-4428-3183

Woojin Kwon https://orcid.org/0000-0003-4022-4132

Hyunju Yoo https://orcid.org/0000-0002-8578-1728

## References

Aspin, C., Reipurth, B., Beck, T. L., et al. 2009, ApJL, 692, L67

Aspin, C., Reipurth, B., Herczeg, G. J., & Capak, P. 2010, ApJL, 719, L50

Bae, J., Hartmann, L., Zhu, Z., & Nelson, R. P. 2014, ApJ, 795, 61

Berry, D. S., Reinhold, K., Jenness, T., & Economou, F. 2007, in ASP Conf. Ser. 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell (San Francisco, CA: ASP), 425

Caratti o Garatti, A., Garcia Lopez, R., Scholz, A., et al. 2011, A&A, 526, L1

Chibueze, J. O., Omodaka, T., Handa, T., et al. 2014, ApJ, 784, 114

Contreras, Y., Sanhueza, P., Jackson, J. M., et al. 2018, ApJ, 861, 14

Covey, K. R., Hillenbrand, L. A., Miller, A. A., et al. 2011, AJ, 141, 40

Currie, M. J., Berry, D. S., Jenness, T., et al. 2014, in ASP Conf. Ser. 485, Astronomical Data Analysis Software and Systems XXIII, ed. N. Manset & P. Forshay (San Francisco, CA: ASP), 391

Dempsey, J. T., Friberg, P., Jenness, T., et al. 2013, MNRAS, 430, 2534

Dunham, M. M., Evans, N. J., II, Terebey, S., Dullemond, C. P., & Young, C. H. 2010, ApJ, 710, 470

Dunham, M. M., & Vorobyov, E. I. 2012, ApJ, 747, 52

Eden, D., Liu, Tie, Kim, K.-T., et al. 2019, MNRAS, 485, 2895

Evans, N. J., II, Dunham, M. M., Jørgensen, J. K., et al. 2009, ApJS, 181, 321

Feng, S., Beuther, H., Zhang, Q., et al. 2016, ApJ, 828, 100

Fischer, W. J., Megeath, S. T., Tobin, J. J., et al. 2012, ApJ, 756, 99

Gibb, A. G., Jenness, T., & Economou, F. 2013, PICARD: A Pipeline for Combining and Analyzing Reduced Data, Starlink User Note 265, (Hilo, HI: Joint Astronomy Centre)

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

Hartmann, L., & Kenyon, S. J. 1996, ARA&A, 34, 207

Herbig, G. H. 1977, ApJ, 217, 693

Herbig, G. H. 2008, AJ, 135, 637

Herczeg, G. J., Johnstone, D., Mairs, S., et al. 2017, ApJ, 849, 43

Hodapp, K. W., Chini, R., Watermann, R., & Lemke, R. 2012, ApJ, 744, 56

Holland, W. S., Bintley, D., Chapin, E. L., et al. 2013, MNRAS, 430, 2513

Hunter, T. R., Brogan, C. L., MacLeod, G., et al. 2017, ApJL, 837, L29

Johnstone, D., Hendricks, B., Herczeg, G. J., & Bruderer, S. 2013, ApJ, 765, 133

Johnstone, D., Herczeg, G. J., Mairs, S., et al. 2018, ApJ, 854, 31

Kenyon, S. J., Hartmann, L. W., Strom, K. M., & Strom, S. E. 1990, AJ, 99, 869

Kóspál, Á., Ábrahám, P., Prusti, T., et al. 2007, A&A, 470, 211

Liu, T., Kim, K.-T., Juvela, M., et al. 2018a, ApJS, 234, 28

Liu, T., Li, P. S., Juvela, M., et al. 2018b, ApJ, 859, 151

Mairs, S., Bell, G. S., Johnstone, D., et al. 2018, ATel, 11583

Mairs, S., Johnstone, D., Kirk, H., et al. 2015, MNRAS, 454, 2557

Mairs, S., Johnstone, D., Kirk, H., et al. 2017a, ApJ, 849, 107

Mairs, S., Lalchand, B., Bower, G. C., et al. 2019, ApJ, 871, 72

Mairs, S., Lane, J., Johnstone, D., et al. 2017b, ApJ, 843, 55

McKee, C. F., & Offner, S. R. R. 2011, in Proc. IAU Symp. 270, Computational Star Formation, ed. J. Alves et al. (Cambridge: Cambridge Univ. Press), 73

Molinari, S., Pezzuto, S., Cesaroni, R., et al. 2008, A&A, 481, 345

Motte, F., Bontemps, S., & Louvet, F. 2018, ARA&A, 56, 41

Ortiz-León, G. N., Dzib, S. A., Kounkel, M. A., et al. 2017, ApJ, 834, 143

Pillai, T., Kauffmann, J., Zhang, Q., et al. 2019, A&A, 622, A54

Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 594, A28

Reid, M. J., Menten, K. M., Brunthaler, A., et al. 2014, ApJ, 783, 130

Reipurth, B., Aspin, C., & Herbig, G. H. 2012, ApJL, 748, L5

Safron, E. J., Fischer, W. J., Megeath, S. T., et al. 2015, ApJL, 800, L5

Sandell, G. 1994, MNRAS, 271, 75

Schmeja, S., & Klessen, R. S. 2004, A&A, 419, 405

Schuller, F., Menten, K. M., Contreras, Y., et al. 2009, A&A, 504, 415

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

Stutzki, J., & Guesten, R. 1990, ApJ, 356, 513

Svoboda, B. E., Shirley, Y. L., Battersby, C., et al. 2016, ApJ, 822, 59

Tan, J. C., Kong, S., Zhang, Y., et al. 2016, ApJL, 821, L3

Terebey, S., Shu, F. H., & Cassen, P. 1984, ApJ, 286, 529

Traficante, A., Fuller, G. A., Billot, N., et al. 2017, MNRAS, 470, 3882

Urquhart, J. S., Figura, C., Wyrowski, F., et al. 2019, MNRAS, 484, 4444

Urquhart, J. S., König, C., Giannetti, A., et al. 2018, MNRAS, 473, 1059

Urquhart, J. S., Moore, T. J. T., Csengeri, T., et al. 2014, MNRAS, 443, 1555

Vorobyov, E. I., & Basu, S. 2015, ApJ, 805, 115

Whitworth, A. P., & Ward-Thompson, D. 2001, ApJ, 547, 317

Williams, J. P., de Geus, E. J., & Blitz, L. 1994, ApJ, 428, 693

Yoo, H., Lee, J.-E., Mairs, S., et al. 2017, ApJ, 849, 69