Tracking Protostellar Variability in Massive Protoclusters with ALMA: I. Insights from QUARKS and MaMMOtH

TL;DR

This study presents the first systematic, multi-epoch exploration of 1.3 mm continuum variability in massive protoclusters using ALMA data from the QUARKS and MaMMOtH surveys, covering 22 protoclusters with timescales from hours to >2 years. A dedicated pipeline aligns images, isolates line-free continuum, and applies relative flux calibration to measure peak-intensity changes of 383 condensations at ~0.3$\arcsec$ resolution. Two complementary variability diagnostics—standard-deviation analysis and difference maps—identify five variable condensations (lower-limit variability fraction ~1.3%), with condensation 10 in I13111-6228 showing a robust ~68% rise over 1 year, consistent with burst-mode accretion scenarios in massive star formation. The work demonstrates the value of high-resolution, multi-epoch interferometry for probing episodic accretion in distant, crowded protoclusters and outlines a path toward larger, longer-baseline variability studies. It also emphasizes limitations due to sampling and beam dilution, motivating future surveys that leverage ALMAGAL overlaps and higher-resolution follow-ups.

Abstract

Millimeter/submillimeter variability is often attributed to dynamical disk-mediated accretion, yet detection is limited to low-mass protostars in nearby clouds. Recent observations have also revealed significant (sub)millimeter variability in high-mass protostars, but the confirmed cases are scarce and lack systematic monitoring. In this work, we analyzed multi-epoch Atacama Large Millimeter/submillimeter Array (ALMA) Band 6 (1.3 mm) continuum observations of 22 massive protoclusters, with epoch separations ranging from a few hours to more than two years, while achieving a consistent angular resolution of approximately 0.3 arcsec. These data allow us to track variability of protostars across a broader mass range and in an environment markedly different from nearby clouds. Using a custom processing pipeline for data reduction, image alignment, and relative flux calibration, we achieve high-precision flux measurements and, for the first time, investigate millimeter variability in massive protoclusters based on interferometric data in a statistical manner. Applying the astrodendro algorithm, we identified 383 condensations and tracked their variations in peak intensities. Standard deviation analysis and difference maps reveal five variable sources, corresponding to a lower limit of 1.3% on the variable fraction. Among these, I13111-6228 stands out as it hosts a hypercompact H II region that exhibits a 68% increase in continuum peak intensity over one year, with an uncertainty of 2%. This supports the burst-mode accretion picture in massive star formation as a viable route for the formation of massive stars.

Figures (15)

• Figure 1: Processing workflow. Schematic overview of the four main processing steps adopted in this study. The process begins with data selection and calibration, including dataset cross-matching and ALMA pipeline calibration to prepare the data for further steps (Sect. \ref{['sec:obs']}). In the second step, imaging and self-calibration, we perform spectral Line flagging and self-calibration (Sect. \ref{['sub:data_flag_selfcalibration']}), followed by smoothing all epochs of each source to a common beam to produce consistent continuum images across epochs (Sect. \ref{['sub:continuum_imaging']}). In the third step, condensation extraction, involves aligning images across epochs, identifying compact condensations in both datasets, and applying a unified mask for consistent source extraction (Sect. \ref{['sub:core_extraction']}). Finally, in flux comparison, we perform relative calibration (Sect. \ref{['sub:relative_calibration']}), track the peak intensity evolution of each condensation across epochs (Sect. \ref{['sub:variables_analysis']}), and evaluate variability using both the standard deviation method (Sect. \ref{['sub:SD_analysis']}) and the difference map analysis (Sect. \ref{['sub:difference_map']}).
• Figure 2: ALMA Band 6 continuum images of I15520--5234 obtained from four epochs spanning over two years (with the observing dates indicated in the figure). The cyan dashed circle in each panel marks the 0.5 primary beam FWHM region of the MaMMOtH survey ($\sim$13.24$\arcsec$), which is also used as the reference area for image alignment. Upper panels (a, b): MaMMOtH 1.3 mm continuum emission. Lower panels (c, d): QUARKS 1.3 mm continuum emission. For consistency, all contour levels in the four panels are defined using the rms noise level measured from the first-epoch image (2022 May 30), where $\sigma_{\rm rms} = 1.23~\mathrm{mJy~beam^{-1}}$. White contours show [5, 10, 20]$\times\sigma_{\rm rms}$, while gray contours show [25, 50, 80]$\times\sigma_{\rm rms}$. The synthesized beam of each epoch is shown in the lower-left corner of each panel, and a 0.05 pc scale bar is indicated in the lower-right corner.
• Figure 3: Top: Cartoon illustrating the extraction of structures and their fluxes. Red and green contours represent condensations independently identified in the MaMMOtH and QUARKS surveys using the astrodendro algorithm. The union of these two sets defines a common mask, shown in yellow, which is consistently applied across all observing epochs. Bottom: Application to I15520--5234. (a) 1.3 mm continuum image before primary beam correction. The dashed black contours indicate levels of [10, 24, 48, 96]$\times\sigma_{\rm rms}$, where $\sigma_{\rm rms}=0.58~\mathrm{mJy~beam^{-1}}$. Red and green contours trace condensations identified independently in the MaMMOtH and QUARKS datasets. (b–e) Continuum images after primary beam correction from 2022 May 30, 2022 August 09, 2024 June 02, and 2024 June 03. The yellow contour traces the union mask used to extract peak intensities for the flux comparison analysis in Section \ref{['sub:SD_analysis']}.
• Figure 4: Relative flux calibration overview. Top: distribution of flux calibration factors. Bottom: corresponding uncertainties.
• Figure 5: Normalized standard deviation of peak intensity ($\mathrm{SD}/\mathrm{SD}_{\mathrm{fid}}$) as a function of the mean peak intensity for condensations across 22 protoclusters. Horizontal dashed lines indicate levels at 1, 3, and 5 times the fiducial expectation ($\mathrm{SD}/\mathrm{SD}_{\mathrm{fid}} = 1, 3, 5$). Sources exceeding the 5 times fiducial level are highlighted in red and flagged as candidate variables.