Young Faithful: The Eruptions of EC 53 as It Cycles through Filling and Draining the Inner Disk

SCUBA-2 (Holland et al. 2013) continuum observations at 450 and 850 μm were obtained simultaneously as part of the JCMT Transient Survey (project code: M16AL001; Herczeg et al. 2017). The Pong 1800 observing mode (Kackley et al. 2010) was used to yield a map with a uniform sensitivity over a 30′ diameter circular region centered at (R.A., decl.) = (18:29:49, +01:15:20, J2000). This region was observed approximately monthly (depending on weather conditions and source availability) since the first observation on 2016 February 2 (UTC). The integration time of each observation was varied to ensure a consistent background noise level of ∼10 mJy/beam at 850 μm (Mairs et al. 2017). The 450 μm background noise level varied greatly as it is much more dependent on submillimeter emission from the Earth’s atmosphere than the 850 μm data. The 450 μm and 850 μm beam sizes were 96 and 141, respectively (Dempsey et al. 2013).

The JCMT submillimeter flux measurements were converted to magnitudes using the following scaling:

where Z₈₅₀ = 1.04 Jy/beam and Z₄₅₀ = 3.58 Jy/beam.

In the near-IR, JHK, EC 53 appears as a compact but resolved source (with an FWHM of around 2″ in the H band) and a nebula that fans out southeastward (e.g., Hodapp 1999). The compact source must contain both the variable EC 53A and its nearby (03 distant) apparently constant companion EC 53B (Hodapp et al. 2012). The nebula can be traced out to about 20″ in 360 s of H-band exposure with the Liverpool Telescope, and presumably extends all the way to the terminal shock S11 an arcminute away (Herbst et al. 1997; Hodapp et al. 2012).

We observed EC 53 using the 3.8 m UKIRT, Hawaii, and the UKIRT Wide Field Camera (WFCAM; Casali et al. 2007). WFCAM employs four 2048 × 2048 HgCdTe Hawaii-II arrays, each with a field of view of 1365 × 1365 at an image scale of 04 per pixel. The target was placed on array #3 of WFCAM during all observations. Monitoring observations were obtained at the J, H, and K bands during a period of over 5 yr from 2014 October 12 to 2020 February 28.

We also obtained J-, H-, and K_s-band imaging photometry of EC 53 with the IRIS 0.8 m telescope (Hodapp et al. 2010) and 1024 × 1024 2.5 μm IRIS IR camera of Universitätssternwarte Bochum on Cerro Armazones, Chile. The observations were part of a long-term monitoring program of the Serpens NW star-forming region. In the J-, H-, and K_s-band filters, the integration time per exposure was 20 s, and for each observation, 10 such frames were obtained, as well as 10 frames on a separate sky position. The data reduction pipeline produced dark-subtracted, flat-fielded, sky-subtracted, and co-added images with a pixel scale of 0375 per pixel. The photometry was extracted using the IRAF task APPHOT with an aperture of 22 pixels = 825 diameter. The photometry was calibrated against a set of Two Micron All Sky Survey stars selected to be isolated, sufficiently bright, and non-variable.

We obtained further IR observations of EC 53 using IO:I (Barnsley et al. 2016) on the Liverpool Telescope (Steele et al. 2004) on the island of La Palma in the Canary Islands. The Liverpool Telescope observed a full light curve for one periodic cycle. Given the extended source structure it is not clear what size aperture is optimal for IR photometry with the Liverpool Telescope. We therefore experimented with 13, 37, and 52 radius apertures using the H-band data. By subtracting the magnitudes in these different apertures we found that the variations between them were less than ±0.05 mag rms. We therefore settled on using a 37 radius aperture, which is greater than the worst seeing but still minimizes the effect of the nebula, for all our IR data. To bring the heterogeneous data sets into alignment, we added an offset to all Liverpool H-band data points. The offset value was optimized by the minimum string length method (Dworetsky 1983) with the combined light curve of UKIRT H and Liverpool H after adding the offset (see Section 3.1 for details about the minimum string length method).

WISE (Wright et al. 2010) is a 40 cm telescope in a low Earth orbit that has surveyed the entire sky using four IR bands centered at 3.4, 4.6, 12, and 22 μm (bands W1, W2, W3, and W4) and with an angular resolution of 61, 64, 65, and 120, respectively. The orbit of WISE allows it to cover every part of the sky at least eight times (Mainzer et al. 2011), with each patch of sky observed many times over a period of ∼1 day. The original survey ran between 2010 January and September. Once the telescope’s cryogen tanks were depleted it continued to operate for 4 months using the W1 and W2 bands (known as the NEOWISE Post-cryogenic Mission; Mainzer et al. 2011). The NEOWISE mission was reactivated in 2013 (Mainzer et al. 2014) and has continued to operate since.

In this work we explored the NEOWISE single-exposure database (2019 data release), which contains W1 and W2 observations from 2013 December to 2018 December (Cutri et al. 2015). The NEOWISE single-exposure detections are complete up to W1 = 15 mag and W2 = 13 mag (Cutri et al. 2015).

To quantify the time dependence of the observed brightness variations, we calculate the period of each light curve using the phase-folded diagram string length method. The string length method (Dworetsky 1983) searches for the period that produces the smoothest phase-folded diagram waveform by minimizing the string length, where the measured length is the sum over the line segments connecting successive points in phase order. We modify the string length method to search for the best-fit offset as well as the best-fit period, recognizing the long-term brightening of the underlying system of EC 53.

We compensate for the overall brightening of the object by increasing the brightness obtained earlier than 2016 September 9 and after 2019 October 4. For each near-IR band, we test brightness offsets between 0 and 1.5 mag, with 0.01 mag intervals, and periods from 500 to 600 days, with 5 day intervals. The modified string length method yields a best fit of 530 days each for J, H, and K. We find 570 days for the 850 μm observations, which are likely affected by the time baseline (see Section 3.4 and Appendix B).

For the data points before 2016 September 9, the best-fit offsets are 1.01, 0.96, and 0.78 mag for J, H, and K, respectively. For the data points after 2019 October 4, we optimize the offset value required to bring the phase curve into agreement. The applied offsets are 0.6, 0.55, and 0.44 mag for J, H, and K, respectively. At 850 μm the offsets are 0.74 mag before 2016 September 9 and 0.24 mag after 2019 October 4, which are applied after the scaling of 5.5 shown in Figure 1.

For completeness, in Appendix B we present a more thorough analysis of the light curves by using periodogram and autocorrelation methods, with and without brightness offsets for the data obtained prior to 2016 September 9 and after 2019 October 4.

In this section we use only epochs with corresponding JHK measurements in order to evaluate the color of EC 53 as a function of phase. Following our results in Section 3.1, we draw the phase-folded diagrams of the light curves in the near-IR and at 850 μm with a 530 day period. We also present the near-IR color as a function of phase. We adopt the same period at every wavelength because the periodicity of the system should be independent of wavelength. Given that our analysis covers roughly three full periods, an uncertainty of less than a few percent in the phase ( days) introduces little phase scatter. The phases are aligned such that the observed minimum in the near-IR occurs at phase 0.5. Our goal in this section is to uncover gross changes in the physical properties of the system with phase.

3.2.1. Phased Light Curves

In the near-IR, the brightness sharply increases starting at phase 0.5, reaches the peak at phase ≈0.6, and decays slowly (Figure 2; note that in the figure we show two complete periods). In the submillimeter band, the brightness increase and decrease appear much more symmetric. Although the submillimeter data points have a large scatter, the initial submillimeter decay stops earlier (phase ≈0.25) than the near-IR drop. The submillimeter rise occurs at roughly the same time as the near-IR rise (phase ≈0.5), but more gently than the near-IR rise; this might suggest that the submillimeter peak is delayed with respect to the near-IR.

Zoom In Zoom Out Reset image size

A significant phase-dependent color variation is also observed (Figure 3). A reddening of the source begins at the same time that the submillimeter curve stops declining (phase ≈0.25), leading to a sharp reddening peak and decline around phase ≈0.6, during the observed brightening of EC 53 in the near-IR. After this abrupt reddening event the color of the source gradually becomes bluer until phase ≈1.

Zoom In Zoom Out Reset image size

The well-sampled color–phase curve reveals an extreme reddening coincident with the abrupt start of each burst. Independently of phase, the near-IR colors of EXors and FUors span a wide range, with J − H = 0.5–3. The color variation observed for EC 53 is distinct from the typical trend found in EXors, which have bluer near-IR colors around the burst and redder colors during the quiescent phase (Kóspál et al. 2011; Lorenzetti et al. 2012). For FUors the observed color variations are more diverse. HBC 722 shows color variation similar to that of EXors (Kóspál et al. 2011) even though the source is classified as an FUor in recent photometric (Connelley & Reipurth 2018) and spectroscopic studies (Lee et al. 2015). The colors of V1057 Cyg and V960 Mon follow an extinction law in the optical band during the decay, right after the peak brightness (Kopatskaya et al. 2013; Hackstein et al. 2015). On the other hand, V346 Nor becomes redder during its slow brightening (Kóspál et al. 2020).

3.2.2. Phased Color–Magnitude Diagrams

To diagnose the cause of the flux and color variation we produce color–magnitude diagrams for the near-IR (Figure 4). Starting at phase ≈−0.25 (green points), the source becomes sharply fainter and slightly bluer. From phase ≈0.25 to phase ≈0.5 (purple and yellow points) it continues to fade but becomes much redder, roughly following the reddening law in typical interstellar matter (Cardelli et al. 1989). The reddening trend in the H band (lower panel of Figure 4) appears flatter than the J-band photometry. Such a flattened trend occurs when the source is reddening and intrinsically brightening simultaneously, which gives insight into the timing offset of the brightening between the near-IR and the submillimeter band (discussed in Section 4). Around phase ≈0.5 (yellow and red points) the source rapidly increases in brightness while staying red. After the abrupt rise in brightness, the source continues to brighten and get bluer, following the reddening law in reverse (red and green points). The source then begins to dim once more (blue points) as the cycle repeats.

Zoom In Zoom Out Reset image size

To measure the decay and rise timescales, we fit a linear function to the decaying and rising sections of the Liverpool Telescope H-band light curve (Figure 5). The measured slopes are 1.47 ± 0.01 mag yr⁻¹ for the dimming and −11.9 ± 0.1 mag yr⁻¹ for the brightening. We convert the slope, n, to the timescale of luminosity variation, τ, as follows:

Here, m₀ and L₀ are the initial magnitude and luminosity, respectively. Thus,

and to determine the timescale in which the luminosity changes by a factor of e,

This reduces to the convenient form, τ ≈ 1.086/n. Therefore, the decay timescale of EC 53 in the H band is ≈0.74 yr, an order of magnitude longer than the brightening timescale of ≈0.091 yr.

Zoom In Zoom Out Reset image size

We overlay the H-band rise and decay slopes measured for a single cycle with the Liverpool Telescope on the full H-band light curve in Figure 6. Each individual rise and decay is well fit with the value derived above, but the timing between the rise and decay shows some variability, indicating an inherent fluctuation within the apparent periodicity of the system.

Zoom In Zoom Out Reset image size

As discussed in Appendix B, the period of EC 53 does not clearly converge to a single value but depends on the time when the observations are taken. The phase-folded diagrams (Figures 2 and 3) also show apparent scatter with respect to a single period.

Measured from Figure 6, the offset times between the observed fadings of EC 53 are 450, 760, and 550 days (blue lines), while the rises are separated by 460, 560, and 550 days (red lines). The time displacement between the crossing points (black squares in Figure 6) indicates the hypothetical maximum and minimum points, with displacements of 459, 582, and 550 days (left to right).

The variations in the mid-IR and 450 μm light curves are consistent with those seen at 850 μm (Figure 7, left panel), though with fewer data points. We perform an orthogonal distance regression to derive the scale factor and offset that give that best agreement between the magnitudes of two different wavelengths. Between the 850 μm and W2 bands, we uncover a scale factor of 6.1 ± 0.5 and an offset of 8.44 ± 0.01 (Figure 7, upper right panel). The best fit for the 850 μm and 450 μm magnitudes yields a scale factor of 2.2 ± 0.1 and an offset of −0.2 ± 0.004 (Figure 7, lower right panel).

Zoom In Zoom Out Reset image size

Contreras Peña et al. (2020) measured a scaling factor of 5.5 ± 0.3 between the submillimeter and mid-IR magnitudes in a sample of a dozen protostars, including EC 53, which simultaneously varied at both wavelengths. The fitting from a single burst of EC 53 was 6.18 ± 0.65. In that paper, the submillimeter to mid-IR correlation was interpreted as the variation of mid-IR emission roughly tracing the underlying bolometric luminosity of the source, and the varying submillimeter emission tracing more closely the dust temperature within the envelope. Additional evidence in favor of the luminosity/dust temperature explanation for the tight correlation and scaling has been found through radiative transfer analyses and spectral energy determinations of modeled embedded sources, both analytic and numeric, with varying central luminosities applied (MacFarlane et al. 2019a, 2019b; Baek et al. 2020).

The measured stronger response to brightness changes at the shorter submillimeter wavelength is also expected, as the 450 μm emission probes somewhat closer to the peak of the dust spectral energy distribution produced by the protostellar envelope (see, for example, MacFarlane et al. 2019a, 2019b; Baek et al. 2020). The detailed radiative transfer modeling in those papers, however, reproduced a scale factor <2, suggesting that there are still aspects of the protostellar models that require refinement.

In Section 3.3 the exponential decay rate of the near-IR brightness is measured as 0.74 yr for the most recent decline, using the Liverpool Telescope data set (Figure 5). This decay rate fits very well the previous decays after burst (Figure 6). Thus, it appears that the timescale has a physical origin, which we associate with a viscous spreading time in the accretion disk surrounding the protostar, . Within order unity (Hartmann et al. 1998),

where

and R₀ is the characteristic size of the accreting disk, H₀ is the scale height of the disk at R₀, c_s is the sound speed in the disk, and is a scaling parameter for the viscosity (see, e.g., Shakura & Sunyaev 1973; Matsuyama et al. 2003).

Rearranging, we obtain an equation for the characteristic size of the accreting disk:

After substitution of fiducial values, we find

Thus, it is possible to associate the decay time with the viscous time in the innermost region of the circumstellar disk, where enhanced accretion is taking place. This requires, however, both a very small inner accretion disk, (10 R_⊙), and a very high viscosity parameter, α ∼ 0.3. We note that this α value is higher than the typical value accepted for FU Ori, ∼0.2–0.02 (Zhu et al. 2007).

The amount of material accreted onto the protostar during a single burst can be determined directly from the bolometric light curve, under the assumption that the observed brightness is due to released accretion energy (Hartmann 2009). Taking a mass of EC 53 of (Lee et al. 2020) and assuming a radius of , we estimate a mass accretion rate of

where η ≤ 1 captures the uncertainty in the fraction of energy radiated away during accretion. It is anticipated that η ∼ 1 (see Hartmann et al. 2011 for a discussion).

The observed exponential decay rate of the near-IR brightness suggests a similar relationship of luminosity versus time,

with (Baek et al. 2020). Substitution into Equation (9) and integration over the length of the observed decay, yr, we find

With η ∼ 1, the maximum accretion rate is , significantly higher than the values typically measured for Class I sources from direct accretion diagnostics (e.g., White & Hillenbrand 2004; Salyk et al. 2013) or from bolometric luminosities (Dunham et al. 2010; Fischer et al. 2017) and even greater than the expected steady-state accretion value for a young, Class 0/I protostar forming over ∼0.5 Myr.

In the preceding sections we have assumed that a viscous inner accretion disk is responsible for the accretion of material onto the protostar. Thus, applying α-disk constraints (e.g., Hartmann et al. 1998), we estimate the required surface density, , in the inner accreting disk:

Rearranging this yields

Furthermore, assuming steady-state α-disk conditions yields and therefore

Thus, the mass of the inner accretion disk is only 2.3 times as large as the total amount of material accreted during a single episodic event suggesting that an outer disk reservoir is required to replenish the inner accretion disk.

While the inner disk is observed to accrete material onto the protostar episodically, we have shown in the preceding section that there must be replenishment from the outer disk in order for the bursts to recur. If we assume that accretion through the outer disk is steady, then

Thus, even our assumed steady-state mass accretion rate in the outer disk reflects better anticipated Class 0 conditions versus Class I expectations.

Furthermore, the outer disk mass is estimated as by Lee et al. (2020), assuming standard conversions between millimeter-continuum emission and mass. If the outer disk has a surface density , then

Starting with a large disk, au (Lee et al. 2020), we extrapolate this outer disk mass to the scale of the inner accretion disk where au, finding . This result is approximately 130 times the above-determined mass of the inner accretion disk (Section 4.3).

We thus conclude that while the outer disk has plenty of mass available near its inner edge, it is unable to smoothly supply the inner disk with material. An accretion blockage occurs somewhere beyond the inner disk, au, where the inwardly accreting material from the outer disk is held up. The blockage is only temporary, however, and after ≈1.5 yr the stored material is released and quickly drains onto and through the inner disk.

An additional consequence of this model is that the effective α value for the inner part of the outer disk is significantly lower than that for the inner disk. Assuming that the quantity does not vary significantly across the inner and outer disk boundary then the ratio of α is inversely proportional to the local surface density scale factor (see Equation (14)), implying that in the outer disk α ∼ 0.002.

From the near-IR light curve we find an exceptionally fast exponential rise time to the peak brightness: yr. Before the near-IR rise, however, we observe enhanced reddening of the source, which we associate with a strong extinction event (Figure 3). The extinction begins around phase 0.25 and becomes extreme near phase 0.5. The full reddening event persists for ≈0.4 yr, although the most extreme () reddening lasts for only ≈0.1 yr.

If the material partaking in each individual burst is held up in the outer disk and then released, rapidly falling to the inner disk, then the fastest timescale available would be the dynamical time. For a protostar, a 0.1 yr infall timescale corresponds to a distance of R ∼ 0.4 au. This scenario requires both a mechanism for holding back the outer disk accretion at ∼0.4 au, and a release mechanism that streams the material directly into the inner disk. Despite these complexities, an intriguing aspect of this scenario is that the burst periodicity, yr, is only a few times the local Keplerian orbital time in the disk where the material is required to be held up (the freefall time is 0.18× the orbital time).

Alternatively, the viscous and accreting outer disk might be efficient to well inside 0.4 au but might become hung up closer to the measured size of the inner disk, 0.05 au, with the effective viscosity at this boundary decreasing steeply such that α ≪1. In this scenario the increased mass buildup near the inner disk boundary takes place over many local Keplerian revolutions, eventually reaching a tipping point during which the effective viscosity rises sharply such that α ∼ 0.3. We note that the physical plausibility of this scenario requires theoretical validation as the effective viscosity has to change dramatically over very short timescales.

The physical cause for the cycles of bursts and blockages is uncertain. Explanations for accretion bursts include magnetospheric instabilities in the inner disk (e.g., D’Angelo & Spruit 2010), gravitational instabilities near the star (e.g., Bae et al. 2014), and interactions between the disk and a close companion (e.g., Bonnell & Bastien 1992; Nayakshin & Lodato 2012; Muñoz & Lai 2016).

The disk fragmentation scenario (e.g., Vorobyov & Basu 2010) is especially intriguing for EC 53 given the estimated disk mass to protostar mass ratio of ∼0.25. This scenario, however, places the fundamental instabilities in the outer disk and would probably require longer timescales than measured for EC 53, although the migration of clumped material from the outer disk might be responsible for triggering an instability within the inner disk (Nayakshin & Lodato 2012).

Of all the possible explanations, the planet or companion hypothesis suggested by Hodapp et al. (2012) and Yoo et al. (2017) is the only one that would specifically invoke quasi-periodic behavior.

Several other outbursts found in the literature have been analyzed for their triggering mechanisms. The mid-IR outburst of Gaia 17bpi preceded the optical outburst, implying an outside–in event where the disk heated prior to accretion onto the central star (Hillenbrand et al. 2018). The outburst event of the Class I object Gaia 19ajj showed extinction from color variation in the mid-IR (Hillenbrand et al. 2019). FU Ori itself may have a constant replenishment of its inner disk (Liu et al. 2019), leading to a smooth light curve for almost a century. Dwarf novae, which undergo repetitive observable outbursts, allow for monitoring and provide considerable references (Osaki 1996). For example, theoretical models have been used to investigate the hysteresis observed in the time-dependent color–magnitude diagrams of some dwarf novae, attempting to fit simultaneously emission from the star, the accretion hot spot, and the disk with moderate success (Hameury et al. 2020). Unlike the hysteresis seen in the color–magnitude diagram for EC 53 (Figure 4), where the change in color with phase appears to be due to a changing extinction, for dwarf novae the optical color change is interpreted as either a variation in the temperature of a component or a change in the relative brightness between components.

Previous studies typically focus on the burst, while in this paper we also demonstrate the relevance of the blockage. In the schematic picture described here (Table 2 and Figure 8), this blockage is understood as a change in the disk viscosity from high α in the innermost regions to low α in the outer regions. While our derived absolute value of α for the inner disk is significantly larger than normally assumed, its ratio against the estimated outer disk α indeed satisfies this requirement.

An alternative possibility invokes changes in the mass-loss rate or pressure within the wind to destabilize disk accretion. Lee et al. (2015) demonstrated that wind pressure prevents the refilling of the disk when accretion is strongest for the early phase of the FUor burst of HBC 722.

Whatever the specific physical causes of the bursts and blockages, they must explain the rise and decay cycles, and the generic picture of the inner disk that we have derived from analyzing those cycles.

Utilizing the observational evidence presented in the preceding section (Section 3) we build here a simple toy model for the near-IR and submillimeter light curves of EC 53. We begin by assuming that the accretion-dominated bolometric luminosity light curve is composed of two coupled exponential functions:

where ϕ is the phase (between 0 and 1), is the offset required to match the observed phase-folded diagram (Figure 2), and the coefficients in the exponential functions reproduce a decay time of 0.74 yr () and a rise time of 0.2 yr () when scaled to the 1.5 yr period of the system. We further assume that the unattenuated near-IR luminosity scales with the accretion luminosity such that

where the nonlinear scaling is meant to create a slight color variation with changing source luminosity such that the unattenuated light from the system is somewhat bluer when least luminous.

The submillimeter brightness is dominated by emission from dust in the extended obscuring envelope surrounding EC 53, which has a delayed reaction of up to a couple of months to changes in the accretion luminosity due to the light-crossing time (Johnstone et al. 2013) for an envelope size of 10,000 au (Baek et al. 2020). To simplify the computation for this model, we calculate the submillimeter flux by averaging the dust temperature response over δϕ = 0.2 and delaying it by (0.15 yr). The top panel of Figure 9 shows the resultant magnitude light curve at the H band and in the submillimeter band (after rescaling). Due to the smoothing and delay, the submillimeter light-curve extrema lag those in the unattenuated near-IR light curve at all phases.

Zoom In Zoom Out Reset image size

The observational data also reveal a strong reddening event just after the submillimeter emission begins to rise (Figure 3), which we associate with enhanced extinction. For the toy model we represent the extinction event as a localized event, using the following equation:

where the maximum extinction takes place at and γ = 0.1.

The effect of the extinction event on the observed near-IR magnitudes is thus

where the coefficients in front of the exponent reproduce the expected reddening law due to extinction (Cohen et al. 1981).

The middle panel of Figure 9 shows the resulting light curves. The submillimeter rise is now advanced with respect to the H band, entirely due to the addition of the near-IR extinction event. The submillimeter peak remains somewhat delayed due to the imposed propagation lag through the extended envelope. These features in the model light curves are now in agreement with the observations. Furthermore the H-band rise has become significantly steeper than the submillimeter rise, which occurs at the same time as the extinction event. As compared with the observations, the most significant departure in the light curve produced by the toy model is the enhanced dip at the end of the H-band decay. In the observations the decay appears smoother, suggesting that the functional form for the extinction event is not quite appropriate. The bottom panel of Figure 9 shows the H − K color for the toy model, again reasonably reproducing the observations.

The toy model color–magnitude diagrams are reproduced in Figure 10. In both panels we include measurements with a two-week cadence in order to show the time spent in each part of the cycle. The overall appearance of the panels matches the observations, despite significant differences in details.

Zoom In Zoom Out Reset image size

The data reduction was performed using the iterative map-making software makemap (see Chapin et al. 2013 for details), which is part of starlink's (Currie et al. 2014) Submillimetre User Reduction Facility (smurf) package (Jenness et al. 2013). Post-processing was applied to each image in order to align it more precisely with the other images and to apply a relative flux calibration. The alignment and relative flux calibration were carried out by analyzing the positions and fluxes of bright, compact, pointlike objects across the field and making relative adjustments based on measurements of non-variable sources (Mairs et al. 2017). Finally, to mitigate pixel-to-pixel noise, smoothing was performed using Gaussians with FWHMs of 4″ and 6″ for the 450 μm and 850 μm images, respectively. All specifics on the data reduction, image alignment, and flux calibration techniques are described in detail by Mairs et al. (2017) (reduction R3).

The observations were performed by dithering the object to five positions separated by a few arcseconds, and a 2 × 2 microstep for each dither. This resulted in a pixel scale of 02 per pixel for the final mosaics. For all epochs, observations in each filter were obtained with two different individual exposure times per frame: 2 s, 2-coadd and 10 s, 1-coadd in J; 1 s, 2-coadd and 10 s, 1-coadd in H; and 1 s, 2-coadd and 5 s, 2-coadd in K.

The sky conditions were not always clear; many of the observations were obtained in the presence of thin or slightly thicker cloud conditions. Preliminary reduction of the data was performed by the Cambridge Astronomical Survey Unit. The details of the data processing, photometric system, and calibration are given in Irwin et al. (2004), Hewett et al. (2006), and Hodgkin et al. (2009), respectively. Further reduction and aperture photometry were performed using the starlink packages KAPPA and PHOTOM. As the sky conditions were not always photometric, relative photometry of EC 53 was obtained using 40 isolated point sources in the field of the same array in which it was located. Care was taken to avoid comparison stars that exhibited variability during the period of our monitoring. Aperture photometry was done in four different apertures, with radii of 129, 15, 368, and 4″. The magnitudes were tied to the JHK magnitudes of the comparison stars measured in the UKIRT Infrared Deep Sky Survey. Our monitoring spanned a period of four outburst cycles. The light curves with the 368 aperture are presented in Figure 1. The error limits given are the 1σ errors in the mean zero-points of the comparison stars.

Although nominally H-band observations, the red edge of the IO:I H-band system response is defined by the sensitivity cutoff of the Teledyne Hawaii-2RG HgCdTe Array. There is no detailed response for the detector, but it is known that its quantum efficiency is only 50% at 1.72 μm,¹¹ compared with a system throughput of 50% of peak throughput at 1.785 μm for the UKIRT system (Hewett et al. 2006). Given that the IO:I H band has an approximately standard blue cutoff at 1.49 μm, this implies the width of the IO:I H band is roughly three-quarters of that for a standard H band.

Each night of observations consisted of several images (usually nine) taken with the telescope moving 14″ in R.A. and/or decl. between each. These individual images were cleaned using the standard Liverpool Telescope IO:I pipeline described by Barnsley et al. (2016). We measured the positions of all the sources in the field using a mean image from one night, and then carried out aperture photometry at those fixed positions (i.e., “forced photometry”) after allowing for offsets and rotations between images.

Having extracted photometry for each star from each image it is standard practice to use an ensemble of non-variable comparison stars to normalize the flux from each image to allow for differences in seeing (which will change the fraction of light that lies outside the aperture) and transparency. There are two issues that make this problematic for the current data. First the IO:I pipeline removes the night sky emission using a median of the other images of the target taken on that night (a “peers-only” median). While the 14″ throw is sufficient to ensure that a small nebula such as EC 53 is moved far enough to obtain a clean sky measurement, there is a much larger region of nebulosity in the field that suffers poor night sky removal, and therefore poor photometry for those stars. Second, the field is relatively sparse as the extinction is high, and the field of view of IO:I is relatively small (627 by 627), exacerbating further the lack of comparison stars.

We were, therefore, careful to examine the precision of the resulting photometry to characterize these effects. We first used all the measurements for each star to derive a χ² with respect to a constant using uncertainties derived from the scatter between sky pixels (e.g., Naylor 1998). We found that we had to add a systematic uncertainty of 0.03 mag to ensure χ² did not increase with the mean flux of the star. This allowed us to remove stars with a reduced χ² of greater than 10 from the normalization procedure on the assumption they were variables, or suffered from poor background removal. We then combined the measurements of each star on each night to create a nightly mean, and rather than using the formal uncertainties, used the scatter between the measurements to create a standard error. Examining the light curves of the remaining stars, we found that our photometry of EC 53 was limited to a precision of 0.04 mag at H = 12 mag for a nightly mean, with the possibility of systematic changes at this level on a timescale of months.

The ACF is a useful tool to determine the period of light curves regardless of their underlying shape. The ACF can be derived by correlating a set of time-series data with itself, shifted by a regular time interval. We adopt Equation (1) of McQuillan et al. (2013) for calculating the ACF of the light curves. In the case of well-sampled periodic data, a clear periodic pattern appears in the ACFs. The position of the second local maximum in the ACF represents the period of variation. Although the ACF requires a set of regularly sampled data, it provides robust and clear results for periodic signals.

Before computing the ACF, the light curves are regularly resampled to have a constant time interval of 15 days, using a linear interpolation. We calculate the ACF for the light curves at near-IR and 850 μm in three different ways (Figure B1): (1) the entire data set at each wavelength (left panels), (2) a common time baseline with 850 μm (middle panels), and (3) the entire data set at each wavelength after removal of the step function determined by the string length method of Section 3.1 (right panels).

Zoom In Zoom Out Reset image size

We identify the position of the second peak in each ACF as the period of the light curve. The uncertainty of the period is estimated at half of the sampling interval (7.5 days). Under the time window constraint but no offset (Figure B1, middle panel), the ACF periods show reasonable agreement across wavelengths, with periods of ≈550 days at all near-IR bands and ≈570 days in the submillimeter band. Applying the offset from Section 3.1 makes the agreement even better, 540 versus 555 days, respectively (Figure B1, right panel). For the unmodified, full near-IR data set, we find the ACF-determined periods are significantly shorter, 510–525 days (Figure B1, left panel). The results are presented in tabular form in Table 3.

Table 3.  The Result of Period Determination

Bands	LSP	ACF	LSPa	ACFa	LSPb	ACFb
J	,	510 ± 7.5		555 ± 7.5		540 ± 7.5
H	,	525 ± 7.5		555 ± 7.5		540 ± 7.5
K	,	525 ± 7.5		555 ± 7.5		540 ± 7.5
850 μm	⋯	⋯		570 ± 7.5		555 ± 7.5

Notes.

^aData points observed after JD = 2,457,421 are used. ^bData points obtained by UKIRT and the offset given before JD = 2,457,640.

Download table as:  ASCIITypeset image

We also apply an LSP (also called a least-squares spectral analysis; Lomb 1976; Scargle 1989) technique to the observed light curves to determine their periods using the LombScargle task in the timeseries package of astropy (Astropy Collaboration et al. 2013). The LSP is a powerful tool to detect periodic signals in discrete time-series data. The false-alarm probability (FAP) is calculated as described in Baluev (2008). For a detailed description, we refer the reader to VanderPlas & Ivezić (2015) and VanderPlas (2018). The best-fit periods along with the uncertainty in the period, determined through Gaussian fitting to each peak in the frequency domain, are tabulated in Table 3.

The LSP-determined periods are very similar to the ACF periods in all cases. For the full, unmodified near-IR data set an additional long-period signal, ≈4 yr, is found, along with the ≈500 day period (Figure B2, left panel). In a common time baseline with 850 μm but without the offset, this additional long-period component endures (Figure B2, middle panel). Applying the string length determined offset to the data prior to 2016 September 9 (Section 3.1), however, removes the long-term component and brings the determined periods, ≈540 days, into agreement with the submillimeter results, the ACF, and the string length method of Section 3.1 (Figure B2, right panel). We note that the FAPs of the peaks obtained by the LSP are well below the 10⁻⁶% level.

Zoom In Zoom Out Reset image size

B.2.1. Periodogram Window Function

We obtain periodograms of the window function from the data set used for Appendix B.2 following Section 7.3.1 of VanderPlas (2018). As shown in Figure B3, there are notable peaks at ≈1 and ≈3 yr in the window function periodograms of the near-IR bands, reflecting the annual observing window for this source. For 850 μm, only minor peaks at 1 yr, ≈200 days, and ≈15 days appear. The 15 day peak corresponds to our observational cadence. None of the LSP-determined periods for EC 53 appears to be influenced by the window function.

Zoom In Zoom Out Reset image size