The JCMT Transient Survey: Six Year Summary of 450/850 μm Protostellar Variability and Calibration Pipeline Version 2.0

The James Clerk Maxwell Telescope (JCMT) Transient Survey has been monitoring eight Gould Belt low-mass star-forming regions since 2015 December and six somewhat more distant intermediate-mass star-forming regions since 2020 February with the Submillimeter Common User Bolometer Array 2 on board JCMT at 450 and 850 μm and with an approximately monthly cadence. We introduce our pipeline v2 relative calibration procedures for image alignment and flux calibration across epochs, improving on our previous pipeline v1 by decreasing measurement uncertainties and providing additional robustness. These new techniques work at both 850 and 450 μm, where version 1 only allowed investigation of the 850 μm data. Pipeline v2 achieves better than 0.″5 relative image alignment, less than a tenth of the submillimeter beam widths. The version 2 relative flux calibration is found to be 1% at 850 μm and <5% at 450 μm. The improvement in the calibration is demonstrated by comparing the two pipelines over the first 4 yr of the survey and recovering additional robust variables with version 2. Using the full 6 yr of the Gould Belt survey, the number of robust variables increases by 50%, and at 450 μm we identify four robust variables, all of which are also robust at 850 μm. The multiwavelength light curves for these sources are investigated and found to be consistent with the variability being due to dust heating within the envelope in response to accretion luminosity changes from the central source.


Introduction
The advent of sensitive large field-of-view detectors has launched an era of time-domain astronomy with (sub)millimeter single-dish telescopes.These data sets have been used to search for and characterize transient events, such as flares from stars (Mairs et al. 2019;Guns et al. 2021;Naess et al. 2021;Johnstone et al. 2022) and relativistic jets (Fuhrmann et al. 2014;Subroweit et al. 2017;Tetarenko et al. 2017), as well as to monitor variability within nearby Galactic star-forming regions associated with the stellar mass assembly process (Johnstone et al. 2018;Park et al. 2019;Lee et al. 2021).Over the next decade, (sub) millimeter time-domain astronomy is anticipated to play an ever more important role as an essential science mode for FYST (CCAT-Prime Collaboration et al. 2023), AtLAST (Ramasawmy et al. 2022), CMB-S4 (Abazajian et al. 2019), and at slightly shorter wavelengths for potential space-based far-infrared missions (André et al. 2019;Fischer et al. 2019).
Despite the significant advances in monitoring capabilities, calibration of (sub)millimeter observations, especially from the ground, remains challenging.Mairs et al. (2021) analyzed over a decade of James Clerk Maxwell Telescope (JCMT) Submillimeter Common User Bolometer Array 2 (SCUBA-2) continuum imager observations and concluded that the peak flux uncertainty at 850 μm after observatory-based calibrations is 7%, while at 450 μm the uncertainty rises to 17%.
Notwithstanding the dynamic range and sensitivity of the modern large format (sub)millimeter detectors, these calibration uncertainties dominate observations of all but the faintest targets.
To overcome this complication, the JCMT Transient Survey (Herczeg et al. 2017) developed a relative calibration scheme (Mairs et al. 2017a(Mairs et al. , 2017b) ) for 850 μm SCUBA-2 observations that significantly improves, by a factor of 3, the default calibration provided by the observatory.This enhanced automated reduction and relative calibration pipeline has allowed the team to identify many years-long secular protostellar variables within the Gould Belt (upward of 30% of the monitored sample; Johnstone et al. 2018;Lee et al. 2021) to examine closely the protostellar variable EC 53 in Serpens (also known as V371 Ser) undergoing episodic accretion events with an 18 month period (Yoo et al. 2017;Lee et al. 2020;Francis et al. 2022), and to investigate the months-long accretion burst associated with the deeply embedded protostar HOPS 373 in Orion B/NGC 2068 (Yoon et al. 2022).Combined, these variability results obtained in the submillimeter and for optically enshrouded protostars have enhanced our understanding of the stellar mass assembly process (Fischer et al. 2022).
The JCMT Transient Survey began as a 3 yr Large Program monitoring eight nearby Gould Belt star-forming regions in 2015 December and has been continually extended without interruption through at least 2024 January.In 2020 February, six slightly more distant intermediate-mass star-forming regions were added to the monitoring observations.Given the large accumulation of time-domain data by this program, we revisit the alignment and calibration strategy of the JCMT Transient Survey to ensure the best quality measurements.The new calibration procedure introduced here also allows for relative calibration of the SCUBA-2 450 μm data sets.Section 2 introduces the observations and the standard data reduction used to make the observatory-calibrated maps at each epoch.Section 3 presents the new image alignment procedure, while Section 4 details the updated relative flux calibration and includes an independent consistency check on its precision.Section 5 presents the coadded images of the eight Gould Belt regions at 850 and 450 μm.We make these deep images publicly available to the astronomical community as part of this paper.We then present a reanalysis of source variability within the Gould Belt regions using the updated data sets in Section 6 and summarize the paper in Section 7.

Observations and Data Reduction
The JCMT Transient Survey (Herczeg et al. 2017) is a James Clerk Maxwell Telescope Large Program (project codes: M16AL001 and M20AL007) dedicated to monitoring the evolution of mass accretion in galactic star-forming regions.The survey employs SCUBA-2 (Holland et al. 2013), a workhorse continuum instrument operating simultaneously at 450 and 850 μm with beam sizes of 9 6 and 14 1, respectively.Since 2015 December, circular fields of ~¢ 30 usable diameter have been obtained approximately monthly, targeting eight star-forming regions in the Gould Belt: IC 348, NGC 1333, NGC 2024, NGC 2068, OMC 2/3, Ophiuchus Core, Serpens Main, and Serpens South (see Herczeg et al. 2017 for details).Combined, these fields contain more than 300 Class 0/I/Flat and more than 1400 Class II young stellar objects (YSOs) as identified by the Spitzer Space Telescope (Megeath et al. 2012;Stutz et al. 2013;Dunham et al. 2015).All data obtained since the beginning of the survey (2015 December 22) to 2022 March 1 are included in this work (the observations are summarized in Appendix A).
In 2020 February, the survey was expanded to monitor six additional fields toward regions of intermediate-/high-mass star formation: three fields in DR21 (north: W75N, central: DR21(OH), south: DR23; see Schneider et al. 2010), M17, M17 SWex (Povich et al. 2016), and S255 (Chavarría et al. 2008).These fields are at distances less than 2 kpc, at which the JCMT beam at 450 and 850 μm are still capable of resolving sub-parsec dust condensations.Previously, these fields have been found to host variable young stars showing evidence of accretion outbursts (Liu et al. 2018;Park et al. 2019;Chen et al. 2021;Wenner et al. 2022).These fields are representative of intermediate-to high-mass star formation regions, in which the earliest high-mass protostars are still in the making.The submillimeter monitoring observations for these fields should, therefore, shed light on the high-mass star-forming process.As they have only recently been added to the survey, these fields have significantly fewer observations than the original eight fields, requiring a separate series of publications for the evaluation of variability.For completeness, a summary of the observations of these six new fields is included in Appendix B.
The exposure time of each JCMT Transient Survey region observation is adjusted based on the amount of atmospheric precipitable water vapor observed along the line of sight to ensure a consistent background rms noise of ∼12 mJy beam −1 at 850 μm from epoch to epoch.The atmospheric absorption is much more severe in the 450 μm transmission band, and the data are more susceptible to atmospheric variability than their 850 μm counterparts.Typical rms noise measurements at 450 μm, therefore, vary over an order of magnitude, between 100 and 1000 mJy beam −1 , for observations that simultaneously yield an uncertainty of ∼12 mJy beam −1 noise at 850 μm (see Section 4.2).Mairs et al. (2017b) describe in detail the data reduction procedures used to construct individual JCMT Transient Survey images.This work continues to make use of the configuration labeled as R3, focusing on the recovery of compact, peaked sources at the expense of accurate extended-emission recovery on scales larger than 200″.Briefly, the MAKEMAP procedure (Chapin et al. 2013), part of the STARLINK software suite's (Currie et al. 2014) SMURF package, is used to iteratively reduce the raw SCUBA-2 data and construct images.In order to well sample the beam, 450 μm maps are gridded with 2″ pixels, while 850 μm maps are comprised of 3″ pixels.A Gaussian smoothing is then performed on each image using a full width at half-maximum equivalent to the angular size of 2 pixels in order to mitigate pixelto-pixel noise, yielding more reliable peak flux measurements.Mairs et al. (2017b) also developed the first version of our procedures for the relative image alignment and flux calibration required to bring individual epochs of the same region into agreement, denoted pipeline v1, or the point-source method.In this work (pipeline v2), we revisit these tasks in order to produce more accurate and robust relative astrometry and calibration.Table 7 in Appendix A summarizes the dates, scan numbers, and rms map noise at both wavelengths, along with the alignment and calibration parameters for all Gould Belt observations.Also, to enhance the value of these JCMT observations, deep coadds of the eight Gould Belt regions are released along with this paper (Section 5).

Image Alignment
The JCMT has an inherent pointing uncertainty of 2″-4″ (see Mairs et al. 2017b and Figure 1).This pointing offset is the same at both 450 and 850 μm because both focal planes have the same field of view on the sky, and observations are carried out simultaneously by means of a dichroic beam splitter (Chapin et al. 2013).Relative image alignment can, therefore, be achieved by determining positional offsets between epochs using only the 850 μm images, for which the background measurement noise is consistent over time (see Section 2).
As part of the JCMT Transient Survey data reduction process for relative alignment, the original point-source method (Mairs et al. 2017b) tracked a set of bright, compact, peaked sources identified in each epoch and compared their measured peak positions to the first image obtained for that region.For each additional epoch, the calculated average offsets in both R.A. and decl.over this set of sources were used to correct the newly obtained data to a fixed grid (Mairs et al. 2017b).While this method produced maps in relative alignment to within a factor of ∼1″, the technique requires at least several bright, pointlike sources to be present throughout the map, and the derived offsets are subject to Gaussian fitting uncertainties introduced by pixel-topixel noise.Furthermore, this method did not attempt to remove the pointing uncertainty of the first epoch.
A more robust pipeline v2 algorithm has been developed and is now included in the updated data reduction pipeline.This new technique cross correlates the reconstructed images associated with each epoch, prior to Gaussian smoothing, to determine the relative alignment offsets and estimates the absolute pointing from the statistics of the individual pointings, assuming that there is no systematic pointing offset at the telescope.In detail, the new image alignment procedure defines a nominally ¢ ´¢ 20 20 subfield centered on the middle of a given map15 taken at epoch i, Map i , and cross correlates the information contained within against the same subfield extracted from the first observation of the given region, Map 0 .The method employs SCIPYʼs 2D, discrete Fourier transform algorithm, FFT2 (Virtanen et al. 2020): where XC is the computed cross-correlation function,  represents the 2D discrete Fourier transform, and * represents the complex conjugate.Map i and Map 0 contain the same structured emission at slightly different positions, tracing the inherent pointing uncertainty between the epochs.Meanwhile, the pixelated random measurement noise is uncorrelated between epochs, and therefore, does not produce a localized peak in XC.The epoch offset is measured by fitting a Gaussian function to XC with a Levenberg-Marquardt least squares algorithm.The necessary offsets that must be applied to Map i in order to bring it into relative alignment with Map 0 are derived by measuring the misalignment of the peak of the cross-correlation function from the center in both the R.A. and decl.directions.In the final step, the central coordinates of Map i are also uniformly shifted to the median R.A. and decl.measured over the unaligned maps observed on or before 2021 June 15, as a reasonable estimate of the true astrometry, assuming pointing errors over many epochs are unbiased.The original point-source method continues to be used as a consistency check for this more robust pipeline v2 algorithm.
The left panel of Figure 1 shows the original telescope pointing offsets with respect to the derived true astrometry toward Orion OMC 2/3, as measured by the XC method, for all epochs (gray squares), the residual offsets measured after the cross-correlation alignment corrections (pipeline v2) have been applied (blue circles), and the residual point-source method alignment consistency checks after the cross-correlation alignment corrections have been applied (yellow triangles).The right panel of Figure 1 shows a zoomed-in comparison between the residual offsets measured by each method.There is excellent agreement between the new and old techniques.The standard deviation of the residual cross-correlation offsets suggests the maps are self-consistently aligned to better than 0 2. The point-source verification suggests the alignment uncertainty is better than ∼0 5, though the inherent uncertainty using the point-source method is larger than for the crosscorrelation method.Table 1 shows the alignment results for each of the eight Transient Survey Gould Belt regions.The ΔR.A. and Δdecl.values are derived by calculating the standard deviation of the measured map offsets from the median (defined to be 0, 0).The corrected residual offsets in the table refer to the point-source alignment verification (yellow triangles in Figure 1), which provides a more conservative measurement of the corrected pointing uncertainty than the cross-correlation self-consistency check.Derived pointing offsets for each Gould Belt observation are given in Table 7 in Appendix A.

Relative Flux Calibration
Both the 450 and 850 μm images are flux calibrated in a relative sense over time by measuring peak fluxes of bright, compact submillimeter emission sources that were originally cataloged in coadded images of each target field.The FELL- WALKER (Berry 2015) source identification algorithm was used to find the locations of compact, peaked sources; it is part of STARLINKʼs (Currie et al. 2014) CUPID (Berry et al. 2013) software package.Sources below brightness thresholds of 1000 mJy beam −1 at 450 μm and 100 mJy beam −1 at 850 μm are excluded from this final source calibration position catalog.Furthermore, sources that were previously identified as known submillimeter variables by Mairs et al. (2017a), Johnstone et al. (2018), and Lee et al. (2021) are removed.
The peak flux of each source in a given catalog is measured in each observed epoch to construct initial light curves for all identified objects.Prior to measuring peak flux, the maps used are Gaussian smoothed by 2 pixels to better match the beam and are relatively aligned to much better than the scale of a single pixel (see Sections 2 and 3).The peak flux in each epoch is then accurately measured at a fixed pixel location without requiring a Gaussian fit for each epoch.
The fiducial (expected) standard deviation, SD fid , in the light curve of a given source, i, has the form (Johnstone et al. 2018 where n rms is the typical rms noise measured across the epochs (dominating faint sources), σ FCF is the expected relative flux calibration uncertainty that can be achieved by the algorithms described below (dominating bright sources), and f m (i) is the mean peak flux of source i.At 850 μm and using the original pipeline v1 calibration approach, Johnstone et al. (2018) found n rms = 12 mJy beam −1 and σ FCF = 0.02 (that is, 2%).

Iteratively Weighted Source Calibration
The pipeline v2 relative flux calibration proceeds as follows.Once the initial light curves for all sources in a given region are constructed, an iterative algorithm is used to weight each source either as a favorable or unfavorable calibrator source based on the robustness of each observed map and the constancy of the source light curve over time.With this information in hand, the full set of calibration sources is used to derive a relative flux calibration factor (relative FCF; R FCF ) for each epoch.The R FCF is a multiplicative constant with which to multiply a given epoch in order to bring the map into agreement with the weighted mean flux brightnesses over all sources, as measured in the coadded image of all data obtained on or before 2021 April 10.The details of the algorithm are as follows: 1. Initial R FCF and their associated uncertainties, σ FCF , are set to be 1% and 2%, respectively, for each epoch and each wavelength.Using these initial estimates, an iterative process begins.2. For each source, excluding known variables, the weighted mean flux, f source ¯is estimated as follows:  2) for 850 μm and 5% of the mean flux at 450 μm), the fiducial value is adopted in order to prevent a runaway effect wherein the brightest, most stable sources obtain outsized weights.4. Source flux thresholds of 1000 mJy beam −1 at 450 μm and 100 mJy beam −1 at 850 μm are then applied to excise sources that are too faint to provide a meaningful contribution to the overall flux calibration.5.The R FCF and formal uncertainty for each individual epoch are then determined by applying a source-weighted linear least squares fit between the peak flux measurements in that observation and their weighted mean values calculated over all observations taken before UTC 2021 April 10.The intercept is fixed at the origin, and the calculated slope and its standard deviation yield the best-fit R FCF and its formal uncertainty, σ FCF,formal .In the equations below, s represents a given source, and e (1.5, 2.2) ( 0.3, 0.5) NGC 2024 (1.5, 1.8) ( 0.8, 0.8) NGC 2068 (1.5, 1.9) represents a given epoch: 6.The R FCF uncertainties (σ FCF ) are estimated by calculating the standard deviation in source fluxes measured within the given epoch normalized to their respective mean fluxes across all epochs.Similarly to Step 3, an empirical analysis yielded a minimum threshold for σ FCF of 70% of the formal R FCF uncertainty to prevent runaway events.Therefore, if σ FCF becomes less than 70% of the formal R FCF uncertainty, σ FCF is set to 0.7 × σ FCF,formal .7. Using the newly calculated relative FCFs, R FCF , and their associated uncertainties, σ FCF , repeat Steps 2-6 until convergence (achieved in at most five iterations).Over time, if a new variable source emerges, it will be automatically downweighted in the calibration and identified through its change in fractional weight as in the left panel of Figure 2 or through the calibration consistency check, described below in Section 4.3.
If all sources have the same fractional uncertainty then they each contribute equal weight.For nonvarying bright sources, this is expected (see Equation (2)) as the calibration dominates the peak flux uncertainty.Alternatively, for faint sources dominated by constant noise, the relative weighting varies directly with the source brightness.Occasionally, sources have additional measurement uncertainties, due to the location on the map, neighbor crowding, or potential variability that has not yet been confirmed (and hence the source has not yet been removed from the list).Due to the increased uncertainty, these sources are automatically downweighted by the algorithm (see point 7 above).All known nonvarying protostars are retained in the pool of calibrators, provided the submillimeter flux is above the thresholds.In this way, deviations in the calibration results over time are used to identify newly varying YSOs and previously undiscovered YSOs, such as deeply embedded protostars that were historically too faint to be included in previous surveys.
Taking the OMC 2/3 region as an example, in the left panel of Figure 2, the upper envelope follows the expected weighting relation, with a few sources lying lower due to higher-thanexpected uncertainties in their fluxes signaling potential variability or difficulty in measuring a reliable peak flux value.The right panel of Figure 2 7 in Appendix A.

Selecting Reliable 450 μm Data
As previously discussed, the 850 μm observations obtained by the JCMT Transient Survey have a consistent background rms noise of ∼12 mJy beam −1 .The 450 μm data, however, are much more sensitive to changes in atmospheric water vapor and airmass, causing more than an order of magnitude of spread in background rms values over the duration of the survey (see Figures 3 and 4).Therefore, in conditions of poor atmospheric transmission, even the brightest sources in 450 μm maps will be overwhelmed by the noise, and no reliable peak flux measurements can be obtained.In order to identify the usable 450 μm maps, in the right panel of Figure 4 we plot a proxy for the atmospheric transmission, the opacity of the atmosphere measured at 225 GHz (τ 225 ) multiplied by the airmass of the observation, against the calculated relative-FCF uncertainty (σ FCF ) for each epoch and identify a box within which the brightest 450 μm peak flux values have a sufficient signal-to-noise ratio (S/N) to return an accurate FCF.The 450 μm map is defined as usable if 1. τ 225 × airmass < 0.12 2. Relative flux calibration uncertainty < 5%

Relative Flux Calibration Consistency Check
An independent relative flux calibration algorithm is employed to verify the consistency of the weighted peak flux calibration scheme described above.The method is based on the original point-source method described by Mairs et al. (2017b), which requires identifying families of bright, reliable (nonvarying) compact sources and constructing light curves (flux versus time).The normalized light curves of the nonvarying family sources trace the inherent JCMT calibration uncertainty from epoch to epoch, i.e., if a given map is brighter than the mean, the average deviation from the mean of the family source peak fluxes will yield a flux correction factor for that epoch.The standard deviation in normalized peak fluxes represents the uncertainty.These results are not used to perform calibration of the data sets.This method is used only as a confirmation of the precision of the weighted peak flux calibration scheme.This verification step does not select reliable 450 μm maps in the same manner as described above.Instead, 450 μm maps are individually selected based on the theoretical R FCF uncertainty that might be achieved given the brightness of the sources present in the map that could constitute a family along with the rms background noise of the map itself.At 450 μm, the target R FCF uncertainty is defined to be 5% to match the weighted flux calibration threshold.The theoretical relative flux calibration uncertainty (σ FCF,theory ) produced by a given family of calibrator sources is calculated by FCF,theory ensemble Here, f is the peak flux of the faintest source considered in mJy beam −1 , N is the number of sources with peak fluxes greater than or equal to f, and rms background noise of a given map in millijansky per beam.Therefore, in order to achieve a relative flux calibration uncertainty better than ∼5%, the rms background noise threshold for whether to consider an individual epoch as reliable is To identify the trustworthy families of sources for each region, the same brightness thresholds for source consideration as in the weighted flux calibration scheme described above, 1000 mJy beam −1 at 450 μm and 100 mJy beam −1 at 850 μm, are initially employed.Unlike in the weighted calibration method, however, not all of these sources will be used to verify the map calibration, they are simply used as a pool of potential calibrator sources.These potential calibrators are arranged in ascending flux order and known variables identified by Mairs et al. (2017a), Johnstone et al. (2018), and Lee et al. (2021) are removed.Beginning with the brightest source and sequentially including fainter sources one by one (excluding known variables), ever larger potential families are defined.For each potential family, Equation (6) is used to determine which maps achieve the rms threshold to produce a theoretical R FCF uncertainty of 5%.At 850 μm, all epochs are included in the analysis.At 450 μm, maps obtained in poor atmospheric transmission conditions are excluded.The set of reliable maps in each region generally matches those identified using the method described in Section 4.2 with a few exceptions (see Appendix A).
At both wavelengths, the potential calibrator groups are then narrowed down by optimizing three key parameters: (1) the brightest set of sources that (2) return the lowest σ FCF using (3) the largest number of maps.The final point applies specifically to 450 μm data.For more detail regarding the trade-off between a higher number of sources and a lower σ FCF , see Section 4.2 of Mairs et al. (2017b).
To test the σ FCF,theory values, the light curves of each potential calibrator group member (in the narrowed-down subset) were constructed using only the maps that met the rms threshold criteria described above for the faintest source in the group.As for the iteratively weighted calibration method, the light-curve fluxes are measured at the known peak pixel location of each source at each aligned epoch.The light curves are then normalized by dividing each flux measurement by the median for that source across the fixed set of observations.
To select the final calibrator family (the group of sources that will be used to relatively flux calibrate the data; see Mairs et al. 2017b), the sources in the narrowed-down potential calibrator group are paired with one another in every possible combination for a given family size from two sources to the number of sources in the group.Figures 5 and 6 summarize the results of the full relative flux calibration procedure.The left panels show the distribution of all σ FCF values derived for the eight Gould Belt regions using the iteratively weighted algorithm described in Section 4.1.The 450 μm data included was observed in good weather conditions and had a σ FCF of 5%, as described in Section 4.2.The inherent telescope flux uncertainty can be estimated from the width of each distribution's main peak.The width of the 850 μm peak is 8%-10%, while the width of the 450 μm peak is 15%-20%, in good agreement with the uncertainties derived via analyses of JCMT calibrator observations (Dempsey et al. 2013;Mairs et al. 2021).
The center and right panels of Figures 5 and 6 show the correlation between the pipeline v2 flux calibration algorithm with the point-source method. 16The 850 μm R FCF values for these two methods agree to better than 3%, while the 450 μm σ FCF values agree to better than 7%.Furthermore, Table 7 in Appendix A allows for the determination of typical conversion uncertainties at 850 μm, σ FCF = 0.01 or 1%, and 450 μm, 0.035 < σ FCF < 0.05 or 3.5%-5%.These values can be used along with Equation (2) when estimating the expected standard deviation of any particular source.

Coadded Images
The JCMT Transient Survey's consistent, repeated observing strategy from 2015 December through 2022 February has led to the deepest 450 and 850 μm maps to date of eight ~¢ 30 -diameter Gould Belt fields.Figures 7 and 8 show each cross-correlation-aligned relative flux-calibrated coadded over the same wavelength-consistent color scales, and Table 2 summarizes the number of epochs included in each image along with the background rms.At 450 μm, only the usable maps, as defined in Section 4.2, were included.
The images released in this work should not be used to analyze large-scale (> ¢ 3 ) flocculent gas and dust structures.The data reduction parameters specifically chosen for this work were optimized to accurately recover compact peak fluxes over a large dynamic range at the expense of reconstructing faint, less dense regions.In fact, while SCUBA-2ʼs combined PONG1800 mapping strategy and data reduction procedure allow for some degree of separation of atmospheric and astronomical signal on scales up to ~¢ 10 , ground-and airbornebased submillimeter bolometer data always suffer from spatial filtering on scales larger than the characteristic detector size, leading to flux loss on the edges of large structures (see Chapin et al. 2013;Mairs et al. 2015Mairs et al. , 2017b;;Mairs et al. 2017a for more details).Coadded images are available on the CANFAR website.17

Results
In this section, for consistent comparisons with previous results, we use the methods outlined by Lee et al. (2021) to redetermine the long-term, secular variables separately at 850 and 450 μm in order to test the updated calibrations.The secular variability is based on both linear and sinusoidal fits to the time-domain observations (see Lomb 1976;Scargle 1989;VanderPlas 2018), with the best-fit false alarm probability (FAP) for interesting sources presented in Table 3.Additional information about specific robust secular variables can be found in the appendix in Lee et al. (2021).Note.Candidate variables are denoted in bold font, and robust variables are denoted in italic font.
We start by comparing the variables recovered by the new 850 μm calibration (pipeline v2) against those recovered by the original scheme (pipeline v1), taking the same 4 yr time window as used by Lee et al. (2021;Section 6.1).We next consider the advantage of a longer time window when searching for secular valuables, utilizing 6 yr of monitoring at 850 μm (Section 6.2).Then, for the first time, we look for evidence of variability at 450 μm, using all 6 yr of monitoring (Section 6.3) and briefly investigate the combined results.

Robustness of 850 μm Variables over 4 yr
Utilizing the original calibration methods (pipeline v1; Mairs et al. 2017b) and the first 4 yr of the JCMT Transient Survey observations through 2020 January, Lee et al. (2021) found 18 secular variables.Based on the best-fit Lomb-Scargle periodogram (Lomb 1976;Scargle 1989; VanderPlas 2018) derived periods P, two of these are periodic (P < 4 yr) and 11 are curved (4 yr < P < 15 yr), while five are best fit as linear-slope sources.In all cases, the FAP was required to be <10 −3 (Figure 9, left panel).These secular variables were all protostars, and no prestellar or disk sources were observed to vary.
Using these methods, every light curve has a best-fit sinusoid, regardless of whether it is a periodic, curved, linear, or even a nonvariable source.The FAP value for nonvariable sources, however, will be high, indicating a poor fit.Furthermore, as these period fits assume an underlying sinusoid, the derived values should be taken as representative of the quantitative nature, rather than specific.For timescales shorter than the observing window, the period is a reasonable measure of episodic timescale, whereas for longer derived periods, the uncertainties are necessarily much larger.Given the <10 yr timescale of these observations, there is insufficient data to distinguish linearly varying sources from sources with sinusoidal variability over long periods since sinusoids can be regarded as a straight line near its zeros when the period is long.Variables with periods longer than 20 yr are, therefore, classified as linear, but also assigned a derived period (see Lee et al. 2021;Park et al. 2021 for further discussions).Further observations over the coming years may refine these fits and allow for differentiation in the future.
With our new calibrations (pipeline v2) and using the same 4 yr of data and FAP threshold (Figure 9, middle panel), we find 22 secular variables; three periodic, 14 curved, and five linear sources.Modifying slightly the classification from Lee et al. (2021), we consider those sources with an FAP <10 −5 to be robust secular variables, and those with a 10 −5 < FAP <10 −3 to be candidate secular variables.With this definition, Lee et al.  (2021) recovered 10 robust and eight candidates, while the new calibrations recovered 13 robust and nine candidates.Robust fit parameters are presented in Table 4.
Comparing the FAP values for two sets of detected secular variables (see Table 3), we find that one source, J034356.5 +320050 in Perseus (also known as IC 348 HH 211), 18 classified as (candidate) periodic by Lee et al. (2021), is rejected using the new calibration.Within the Lee et al. (2021) sample, HH 211 had the shortest period, ∼1 yr, and the second-largest FAP = 10 −3.3 .With the new calibration, the sinusoidal fitting finds the same best period; however, the FAP increases significantly to FAP = 10 −0.5 .Of the other seven candidate sources from Lee et al. (2021), four remain candidate variables using the new calibration, while three are elevated to robust detections.Furthermore, all 10 of the Lee et al. (2021) sources with an FAP <10 −5 remain robust with the new calibration.Finally, we note that five additional sources are classified as candidate secular variables using the new calibration, of which three are unassociated with known YSOs.All three of these sources, however, have FAPs within a factor of 2 of the cutoff threshold.
Comparing our results between both of the calibration methods also provides a test on whether the calibration method affects the variability measurements and best-fit parameters.The left panel of Figure 10 plots the derived light-curve linear slopes for the sources that were classified as variables, candidate and robust, for both calibration methods.The slopes are consistent between the calibrations, as expected.There are significant changes for a few of the fitted periods, however.The right panel of Figure 10 shows the comparison of the best-fit periods for those 10 sources with robust periodogram FAPs using the old calibration.Eight out of 10 sources maintain their secular type (periodic/curved/linear), though these sources often recover somewhat different best-fit periods between the calibrations.The other two sources switch between linear and curved type.Figure 11 plots the old (red) and new (blue) calibrated measurements and best fits and shows that the change to the curved type is subtle, depending specifically on the calibration of the early and late epochs.For the slope comparison, all sources with candidate or robust secular fits are included.For the period comparison, only those sources with robust periodic fits are included, and the dashed horizontal and vertical lines separate the plot into periodic, curved, and linear regions.Two sources are overlapped in the period plot, so each is annotated with a bullseye.Figure 11.Light curves of two sources that changed variable type from linear (red, pipeline v1) to curved (blue, pipeline v2) after the improved calibration. 18In Tables 3 and 4 in Lee et al. (2021), the source names for HH 211 and IC 348 MMS 1 were inadvertently reversed.

Recovered 850 μm Variables at 6 yr
Increasing the observing window from 4-6 yr, from the beginning of the survey through 2022 February, we update the criteria for secular variability accordingly.We categorize the variables by their best-fit periods with periodic (<6 yr), curved (6-20 yr), and linear (>20 yr).We now recover 38 secular variables, more than double the number obtained over 4 yr by Lee et al. (2021).Of these sources, 20 are robust and 18 are candidate variables (Table 3).Of the robust sources, eight are linear, 10 are curved, and two are periodic (Table 5).All robust sources are known protostars except the brightening curved source J182947.6+011553 in Serpens Main, which is discussed further below.
Of the candidate sources, five are linear, six are curved, and seven are periodic.Half of the candidate periodic sources have very short periods, less than a year, and FAPs close to the cutoff value.Another seven of the candidate sources are not known to be associated with YSOs.We suspect that systematic issues, potentially associated with the yearly weather patterns at Maunakea, are responsible for some of these candidates.The right panel of Figure 9 plots the sinusoidal versus linear FAPs for all the monitored JCMT Transient sources over the 6 yr window using the updated calibration, revealing a cluster of sources with FAPs of ∼10 −4 .
Considering only the newly calibrated results, we find that 12 of the 13 robust secular variables after 4 yr remain robust after 6 yr of observations, with only Serpens Main Solar Maximum Mission (SMM) 10 missing the cut (Table 3).The FAP for SMM 10 increases from 10 −5.1 to 10 −4.3 , dropping it to candidate status.Interestingly, SMM 10 was found to have significant stochasticity in its light curve when observed with higher angular resolution using the Atacama Large Millimeter/ submillimeter Array Atacama Compact Array (Francis et al. 2022).Furthermore, of the nine candidate secular variables found after 4 yr, four are raised to robust after 6 yr, while three remain candidates and two are removed from the variable sample.
The majority of the robust sources maintain the same periodic/curved/linear classification over both time windows, although the additional epochs do occasionally change the bestfit period.As an example, Figure 12 plots the old (red) and new (blue) calibrated measurements and best fits for a source that changed from linear type to curved type through the addition of two more years of observation.
As noted above, one robust secular variable is not identified with a previously known protostar, J182947.6+011553 in Serpens Main.Considering near-IR variability, Kaas (1999) cataloged a candidate YSO, source 6 in their Table 5, about 9″ away; however, no object is visible as a localized source in the Wide-field Infrared Survey Explorer (WISE) mid-IR images at this position.Additionally, the source is not in the Spitzer 24 or 4.5 μm data, nor the eHOPS catalog in IRSA, suggesting that it may be a variable background galaxy.Looking further afield, SMM 1 (Casali et al. 1993;Enoch et al. 2009), a robust secular 850 μm variable source (Table 5), lies about 45″ to the southeast of J182947.6+011553.SMM 1 is found to be a secular variable  (Hull et al. 2016) in the direction of J182947.6+011553.Consistent with the scattering surface seen in WISE images, theoretical modeling of the outflow by Liang et al. (2020) suggest an ~¢ 2 lobe length, which would entirely surround J182947.6+011553.Finally, both SMM 1 and J182947.6+011553show rising light curves, with an ∼1.5%-2% increase per year (Table 5).Thus, we hypothesize that this source, despite the significant distance, is being influenced by SMM 1.

Recovered 450 μm Variables at 6 yr
Figure 13 plots the 450 μm sinusoidal versus linear FAPs for all the monitored JCMT Transient sources.There are four robust secular detections, all of which are also robust at 850 μm (see Tables 3 and 6).Furthermore, of the 12 additional candidate variables at 450 μm, four are known robust variables at 850 μm.In Figure 14, we show the 450 and 850 μm light curves for the four robust at 450 μm protostellar sources, along with a scatter plot of the paired submillimeter fluxes across all epochs, revealing the linearity in the response between wavelengths.This is expected if the underlying process is a change in the temperature of the dust in the envelope due to variations in the accretion luminosity of the deeply embedded protostar (Johnstone et al. 2013;Contreras Peña et al. 2020;Francis et al. 2022).
Quantifying the response at 850 μm versus 450 μm, we have calculated the normalized slopes for each of these four sources (see Figure 14, right panels).For the three Orion HOPS secular variables, the normalized slope is ∼1.4,such that the 450 μm brightness has about a 40% larger variation than the 850 μm brightness.Following the same argument as used by Contreras Peña et al. (2020, their Section 6.2), we anticipate that the 850 μm brightness varies as Td 1.5 , while the 450 μm brightness varies as Td 2 , where we assume T d ∼ 20 K in the outer envelope.The larger exponent at 450 μm is due to the fact that the shorter wavelength is less fully on the Rayleigh-Jeans tail of the dust emission, and therefore, has a stronger reaction to temperature change.These two formulae can be combined to yield µ S S 450 850 4 3 , which for small variations can be made linear such that S 450 ∼ 1.33 S 850 .
The Serpens Main secular variable protostar EC 53 (Lee et al. 2020;Francis et al. 2022), also known as V371 Ser, shows a significantly higher normalized slope, ∼1.8.This may indicate additional extended and nonvarying emission at 850 μm, which is biasing the slope higher, or that the mean dust temperature T d used in the above determinations of submillimeter brightness response has been overestimated.The power-law exponents rise at both 850 and 450 μm as the dust temperature decreases, but more strongly at 450 μm, leading to a greater response at 450 μm versus 850 μm.We note, however, that modeling of the EC 53 envelope structure and temperature profile by Francis et al. (2022) to fit simultaneously time-variable observations at interferometric and single-dish angular scales required a somewhat higher outer dust temperature, T d ∼ 25 K. Francis et al. (2022) also struggled to fit well these 450 μm JCMT observations (see their Section 6.2) and suggest that one complication may be the JCMT beam structure at 450 μm.The good fits, however, to the simple Contreras Peña et al. (2020) model for the three HOPS sources suggest that the problem may lie with the simplified modeled structure assumed within the EC 53 envelope.Thus, time-dependent calculations using a more detailed and careful radiative transfer modeling of the envelope, such as those considered by Baek et al. (2020) and including the known outflow cavity and disk, are likely to be required.
Finally, as a check on the utility of the robust secular fits at both 850 and 450 μm, in Figure 15, we present scatter plots at 850 and 450 μm of the reduced χ 2 versus mean flux (using the middle 80% of data points in brightness).The measurements are made both before and after removal of the best-fit secular component for those sources that are robust secular variables at 850 μm over 6 yr.As found by Lee et al. (2021), the robust variables tend to have exceptionally high χ 2 values prior to the removal of the best secular fit.For the 850 μm data sets, even the removal of the secular fit can leave a significant residual, suggesting that for some sources, there exists an additional component beyond the simple smooth evolution in time of the light curves investigated here.

Summary and Conclusions
The JCMT Transient Survey has been monitoring eight Gould Belt regions continuously since 2015 December and six intermediate-mass star-forming regions since 2020 February.
We have reduced and calibrated all of these observations through 2022 February, using pipeline v2, which we describe in this paper.We have also investigated the Gould Belt fields for variable sources, comparing the 4 yr results using pipeline v2 against the analysis by Lee et al. (2021) over the same time range, using pipeline v1.We further show the importance of increasing the monitoring observations time window by comparing the variability results using pipeline v2 over the full 6 yr of individual epochs.Finally, we demonstrate the value of multiwavelength monitoring by performing pipeline v2 relative calibration on the 450 μm data and comparing individual source light curves at both 850 and 450 μm.The main results of this paper are therefore as follows: 1.The pipeline v2 relative alignment via cross correlation achieves an accuracy better than 0 5, which can be compared directly against the 850 and 450 μm beam sizes of 14 1 and 9 6, respectively.While the alignment is performed at 850 μm, the results are directly applicable to the 450 μm maps as the JCMT SCUBA-2 observations are carried out simultaneously by means of a dichroic beam splitter (Section 3). 2. The pipeline v2 relative flux calibration uses information on all sources in each field and an iterative approach to more robustly bring all epochs of a given region into agreement.For most fields, the uncertainty in the 850 μm relative FCF (σ FCF ) is 1%, and at 450 μm, the relative-FCF uncertainty is <5% (Section 4).Furthermore, given that the JCMT Transient Survey observations are optimized for 850 μm observations, we automate the determination of the good 450 μm epochs by applying limits to the quantified metadata (Section 4.2). 3. Along with this paper we make available to the community deep coadds of the eight Gould Belt starforming regions (Section 5).The typical rms for these maps at 850 and 450 μm are 1.5 and 24 mJy bm −1 , respectively (Table 2).We also present the reduction and calibration metadata for all Gould Belt epochs (Appendix A) and all intermediate-mass star-forming regions (Appendix B).
4. We analyze the variable source results at 850 μm using both pipeline v2 and pipeline v1 over 4 yr, uncovering a greater number of robust variables using the updated calibration techniques (Section 6.1).The properties of the recovered variables found by both calibration methods agree in general, though there are subtleties when determining the best-fit periods via periodogram analysis.
We also demonstrate that extending the time baseline to 6 yr increases the number of variables recovered (Section 6.2). 5. Finally, we analyze the variable source results at 450 μm for the first time.Over 6 yr, we recover four robust variables, all of which are also robust at 850 μm.Direct comparison of the light curves at both wavelengths supports the expectation that the variability is driven by changing mass accretion onto the central protostar and the resultant heating of the dust in the natal envelope (Section 6.3).
Coadded images are available on the CANFAR website (see footnote 17).

Figure 1 .
Figure 1.Region: OMC 2/3.Inherent positional offsets as measured by the cross-correlation algorithm (this work) are shown as gray squares.Residual offsets, after applying pointing corrections to each image, are represented by blue circles.Yellow triangles show the results of a consistency check performed on the corrected images using the original (pipeline v1; point source) pointing correction algorithm.
where e represents each epoch, N e represents the number of epochs, and f source,e represents the flux of the given source in epoch e.3.The source uncertainty, σ source , is then estimated by calculating the standard deviation in the source flux across all epochs, weighted by the epoch uncertainty.If the calculated source standard deviation is smaller than the fiducial standard deviation (Equation ( plots the resultant R FCF required for each individual epoch.The upward trend shown in the right panel is seen across all Gould Belt regions and corresponds to systematic changes in the JCMT flux throughput since 2015, which changed the standard FCFs published by the observatory, described in detail by Mairs et al. (2021).These changes were (1) a filter stack replacement in 2016 November and (2) secondary mirror unit maintenance that improved the beam profiles in 2018 June.The flux conversion factors presented in Mairs et al. (2021) are reciprocals of those derived in this work and thus decrease over time, while the R FCF presented here appears to increase.The magnitude of the standard-FCF changes agrees in both of these studies.The derived R FCF values for each observation are given in Table

Figure 2 .
Figure 2. Region: OMC 2/3.Left: the fractional weight each identified pointlike source carries in the weighted, relative flux calibration.Right: calculated relative FCFs as a function of time.The vertical line denotes the date before which the median flux is calculated for each source to employ in the relative normalization (UTC 2021 April 10; see point 2 in Section 4.1).The dashed horizontal lines show the 1σ range of the relative-FCF values.

Figure 4 .
Figure 4. Left: 450 (red) and 850 μm (blue) rms vs. τ 225 × airmass (a proxy for the atmospheric transmission).The vertical black line indicates the τ 225 × airmass threshold, as indicated in the right panel.Right: τ 225 × airmass vs. the weighted calibration factor uncertainty (pipeline v2).The blue data points indicate the usable maps as described in Section 4.2.

Figure 5 .
Figure5.All low-mass regions, 850 μm.Left: the distribution of relative FCFs computed using the weighted algorithm.Center: relative FCFs computed by the original pipeline v1 (point source) as a function of the relative FCF computed using the weighted algorithm (this work).A 1:1 line is overlaid.Right: distribution of the point-source algorithm relative FCFs divided by the weighted algorithm relative FCFs, showing consistency between the two methods to within 2%.
For each pair of sources in each family size, the normalized fluxes are compared to one another, epoch by epoch.The standard deviation of this normalized flux ratio is an indication of how well the sources agree with one another as they track the inherent JCMT flux calibration uncertainty from epoch to epoch.Higher standard deviation values indicate the sources are in less agreement.The calibrator family is selected by optimizing the group size with the minimum standard deviation threshold of each group.Finally, in each epoch, this secondary R FCF is used to verify the main weighted calibration method, taking the average normalized flux value among the sources comprising the optimized family.The normalization is calculated over the same time period as the pipeline v2 data taken prior to 2021 April 10 UTC such that the solutions are fixed for future observations.

Figure 6 .
Figure 6.Same as Figure 5, but for the good 450 μm data (see Section 4.2).

Figure 8 .
Figure 8. Coadded images of IC 348, NGC 1333, Serpens Main, and Serpens South (top to bottom, respectively) at both 450 (left) and 850 μm (right).Observations through 2022 February are included.At 450 μm, only usable data as defined in Section 4.2 are included.

Figure 9 .
Figure9.Scatter plot of sinusoidal (FAP Mod ) and linear (FAP Lin ) at 850 μm.Left: 4 yr of data using the original alignment and reduction pipeline as performed byLee et al. (2021).Middle: 4 yr of data reduced with the new calibration pipeline introduced in this paper.Right: 6 yr of data reduced with the new calibration pipeline.Note that the x-axis extends to much smaller FAP values in the right panel compared with the left and center panels.

Figure 10 .
Figure10.Comparison of best-fit linear slope (left) and period (right) for old and new calibrations at 850 μm.For the slope comparison, all sources with candidate or robust secular fits are included.For the period comparison, only those sources with robust periodic fits are included, and the dashed horizontal and vertical lines separate the plot into periodic, curved, and linear regions.Two sources are overlapped in the period plot, so each is annotated with a bullseye.

Figure 14 .
Figure 14.Left panels: light curves at 850 μm (red dots) and 450 μm (blue dots) for the four robust variables at both wavelengths, showing the tight correlation.Right panels: scatter plots of the brightness at each wavelength normalized at 850 μm by the median measured value and at 450 μm by the intercept value of the best-fit slope at the median 850 μm value.The normalized slopes are provided for each source.

Figure 15 .
Figure 15.Scatter plot of the χ 2 -fit value divided by the number of observations for each source.The left figure is drawn at 850 μm, and the right figure is drawn at 450 μm.Here, we only include the middle 80% data points in brightness and use Equation (2) for the standard deviation of each measurement, with RFCF unc equal to 1% and 5% at 850 and 450 μm, respectively.The robust variables at 850 μm are also annotated.The gray-vertical lines show the change in the χ 2 value after the subtraction of the best-fit secular trend.

Table 7 in
Appendix A includes the derived 450 μm R FCF values for those epochs that satisfy the above conditions.The table also allows for a calculation of the typical 450 μm rms noise in an epoch, n rms = 130 ± 65 mJy beam −1 , which can be used in Equation (2) to estimate the fiducial 450 μm standard deviation of any particular monitored source.

Table 2
Summary of Gould Belt Region Coadded Images Released

Table 4
Physical Properties of Robust Secular Variables at 850 μm over 4 yr

Table 5
Physical Properties of Robust Secular Variables at 850 μm over 6 yr Figure 12.Light curve of HOPS 358, whose variable type changed from linear (red) to curved (blue) after adding two additional years of observation.The dashed line marks the date of the last measurement used by Lee et al. (2021) in their linear fit.
in the mid-IR (Contreras Peña et al. 2020; Park et al. 2021), which has launched a powerful jet and outflow

Table 6
Physical Properties of Robust Secular Variables at 450 μm over 6 yr Figure 13.Scatter plot of sinusoidal (FAP mod ) and linear (FAP lin ) at 450 μm.Sources that are robust variables at 850 μm are marked as circles (periodic), triangles (curved), and squares (linear).

Table 7
Summary of Processed Observations of JCMT Transient Survey Gould Belt Regions The entire table is published only in the electronic edition of the article.The first five lines are shown here for guidance regarding its form and content.
(This table is available in its entirety in machine-readable form.)

Table 9
Summary of JCMT Transient Survey Intermediate-Mass Region Processed Observations (This table is available in machine-readable form.)